text:
      PROGRAM CCPSET
C
C Define the error file, the Fortran unit number, the workstation type,
C and the workstation ID to be used in calls to GKS routines.
C
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)   ! NCGM
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=8,  IWKID=1)   ! X Windows
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=11, IWKID=1)   ! PDF
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=20, IWKID=1)   ! PostScript
C
      PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)

      PARAMETER (K=40,N=40,LRWK=1000,LIWK=1000)
      REAL Z(K,N), RWRK(LRWK)
      INTEGER M, IWRK(LIWK)

      CALL GETDAT (Z, K, M, N)
C
C Open GKS
C
      CALL GOPKS (IERRF, ISZDM)
      CALL GOPWK (IWKID, LUNIT, IWTYPE)
      CALL GACWK (IWKID)
C
C Draw a perimeter around the whole frame to show off SET call
C
      CALL SET (0.,1.,0.,1.,0.,1.,0.,1.,1)
      CALL PERIM(0,0,0,0)
C
C Force Data to be drawn in a rectangle twice as long as it is wide
C and only draw in the lower right hand corner of the screen
C
      CALL SET (0.5,0.98,0.125,0.375,0.,2.,0.,1.,1)
      CALL CPSETI('SET - DO-SET-CALL FLAG',0)
      CALL CPSETR('XC1 - X COORDINATE AT INDEX 1',0.)
      CALL CPSETR('XCM - X COORDINATE AT INDEX M',2.)
      CALL CPSETR('YC1 - Y COORDINATE AT INDEX 1',0.)
      CALL CPSETR('YCN - Y COORDINATE AT INDEX N',1.)
C
C Initialize Conpack
C
      CALL CPRECT(Z,K,M,N,RWRK,LRWK,IWRK,LIWK)
C
C Draw perimeter
C
      CALL CPBACK(Z, RWRK, IWRK)
C
C Draw Contours
C
      CALL CPCLDR(Z,RWRK,IWRK)
C
C Close frame and close GKS
C
      CALL FRAME
      CALL GDAWK (IWKID)
      CALL GCLWK (IWKID)
      CALL GCLKS

      STOP
      END

      SUBROUTINE GETDAT (Z, K, M, N)
      INTEGER I,J,K,M,N
      REAL Z(K,N)

      M=K
      DO 10, I=1,M
         DO 20, J=1,N
            Z(I,J)= 10.E-5*(-16.*REAL(I**2*J) +
     +           34*REAL(I*J**2) - REAL(6*I) + 93.)
 20      CONTINUE
 10   CONTINUE

      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.

!> @brief compf(i,j,k,ismpl) calculates the compressibility of the fluid
!> @param[in] i cell index, direction I0
!> @param[in] j cell index, direction J0
!> @param[in] k cell index, direction K0
!> @param[in] ismpl local sample index
!> @return compf  [1 / Pa]
!> @details
!> compf(i,j,k,ismpl) calculates the compressibility in (in 1/Pa) of
!> salin water, given temperature (t, in c) pressure (p,in pa), and
!> salinity (s, in mol/L) at node(i,j,k)\n\n
!>
!> Sources:\n
!> Driesner & Heinrich, The system H2O-NaCl. Part I: Correlation
!>    formulae for phase relations in temperature-pressure-composition
!>    space from 0 to 1000 C, 0 to 5000 bar, and 0 to 1 XNaCl
!>    Geochimica et Cosmochimica Acta 71 (2007) 4880-4901\n\n
!>
!> Driesner, The system H2O-NaCl. Part II: Correlations for molar
!>    volume, enthalpy, and isobaric heat capacity from 0 to 1000 C, 1
!>    to 5000 bar, and 0 to 1 XNaCl Geochimica et Cosmochimica Acta 71
!>    (2007) 4902-4919\n\n
      double precision function compf(i,j,k,ismpl)
        use arrays, only: temp, pres, tsal
        use mod_flow, only: pa_conv1

        implicit none

        ! Location indices
        integer, intent (in) :: i
        integer, intent (in) :: j
        integer, intent (in) :: k

        ! Sample index
        integer :: ismpl

        ! Molar mass of NaCl [g/mol]
        double precision, parameter :: mmnacl = 58.44277d0

        ! Mass fraction of NaCl in solution [-]
        double precision :: fracnacl

        ! Molar mass of water, H2O [g/mol]
        double precision, parameter :: mmwater = 18.01528d0

        ! Temperature, [degC]
        double precision :: t

        ! Pressure [MPa]
        double precision :: p

        ! Monomials of pressure [bar]
        double precision :: pb, pb2, pb3

        ! Pure water compressibility [1/Pa]
        double precision cw
        double precision, external :: compw

        ! Pure water density [kg/m3]
        double precision rw
        double precision, external :: rhow

        ! Local salinity [mol/kg]
        double precision :: s

        ! Mole fraction of NaCl [-]
        double precision :: xnacl

        ! Temperature at which pure water has the same molar volume
        ! as the solution [K]
        double precision :: tv

        ! Driesner2007: Coefficients

        ! Equation (9)
        double precision :: n1, n10, n11, n12

        ! Equation (10)
        double precision :: n2, n20, n21, n22, n23

        ! Equation (11), (12)
        double precision :: n1x, n2x

        ! Table 4
        double precision :: n30, n300, n301, n302
        double precision :: n31, n310, n311, n312

        ! Deviation function (13)
        double precision :: dt

        ! Local temperature [degC]
        t = temp(i,j,k,ismpl)

        ! Local pressure [MPa]
        p = pres(i,j,k,ismpl)*pa_conv1

        ! Local salinity [mol/L]
        s = tsal(i,j,k,ismpl)

        ! pure water compressibility [1 / Pa]
        cw = compw(p,t)

        if (s<=0.0d0) then

          compf = cw

        else
          ! pure water density [g/L = kg/m3]
          rw = rhow(p,t)

          ! Pressure [bar] and pressure monomials
          pb = p*10.0d0
          pb2 = pb*pb
          pb3 = pb*pb2

          ! Mass fraction of NaCl, mol/L > (g/L) / (g/L) mass fraction
          fracnacl = s*mmnacl/(rw+s*mmnacl)

          ! Mole fraction of NaCl, (mol/L) / (mol/L)
          xnacl = (fracnacl/mmnacl) / (fracnacl/mmnacl+(1-fracnacl)/mmwater)

          ! Driesner2007 parameters, Table 4
          n11 = -54.2958d0 - 45.7623d0*exp(-9.44785d-4*pb)
          n21 = -2.6142d0 - 0.000239092d0*pb
          n22 = 0.0356828d0 + 4.37235d-6*pb + 2.0566d-9*pb2
          n300 = 7.60664d6 / (pb+472.051d0)**2
          n301 = -50.0d0 - 86.1446d0*exp(-6.21128d-4*pb)
          n302 = 294.318d0 * exp(-5.66735d-3*pb)
          n310 = -0.0732761d0 * exp(-2.3772d-3*pb) - 5.2948d-5*pb
          n311 = -47.2747d0 + 24.3653d0*exp(-1.25533D-3*pb)
          n312 = -0.278529 - 0.00081381*pb

          ! Driesner2007: n1 and n2 parameters for liquid NaCl
          n1x = 330.47d0 + 0.942876d0*sqrt(pb) + 0.0817193*pb - &
            2.47556D-8*pb2 + 3.45052D-10*pb3
          n2x = -0.0370751 + 0.00237723*sqrt(pb) + 5.42049D-5*pb + &
            5.84709D-9*pb2 - 5.99373D-13*pb3

          ! Driesner2007: Set xnacl=1 in (9), => n1 = n1x
          n10 = n1x

          ! Driesner2007: Set xnacl=0 in (9) => n1 = 0.0d0
          n12 = -n10 - n11

          ! Driesner2007: Set xnacl=0 in (10) => n2 = 1.0d0
          n20 = 1.0d0 - n21*sqrt(n22)

          ! Driesner2007: Set xnacl=1 in (10) => n2 = n2x
          n23 = n2x - n20 - n21*sqrt(1.0d0+n22)

          ! Driesner2007: Equation (9)
          n1 = n10 + n11*(1.0d0-xnacl) + n12*(1.0d0-xnacl)**2

          ! Driesner2007: Equation (10)
          n2 = n20 + n21*sqrt(xnacl+n22) + n23*xnacl

          ! Driesner2007: deviation function D(T) Equations (14)-(16)
          n30 = n300 * (exp(n301*xnacl) - 1.0d0) + n302 * xnacl
          n31 = n310 * exp(n311*xnacl) + n312 * xnacl
          dt = n30 * exp(n31 * t)

          ! Temperature at which pure water has the same molar volume
          ! as the solution
          tv = n1 + n2*t+dt

          ! Pure water compressibility at tv
          compf = compw(p,tv)

        end if

        return

      end function compf


text:
program print_test
    use, intrinsic :: iso_fortran_env, only : SP => REAL32, DP => REAL64
    implicit none
    integer, parameter :: n = 15
    real(kind=SP), parameter :: x = 3.14_SP
    real(kind=DP), parameter :: y = 2.69_DP
    character(len=10), parameter :: c = 'abc'
    logical, parameter :: l = .true.
    type udf
        real :: x
        integer :: n
    end type udf
    type(udf) :: u
    real, dimension(3) :: v = [ 3.5, 5.7, 7.3 ]
    integer, dimension(2, 3) :: A = reshape( [ 3, 7, 5, 9, 2, 4 ], [ 2, 3 ])
    u = udf(x = 3.5, n = 7)

    print *, '|', n, '|'
    print *, '|', x, '|'
    print *, '|', y, '|'
    print *, '|', c, '|'
    print *, '|', l, '|'
    print *, '|', u, '|'
    print *, '|', v, '|'
    print *, '|', A, '|'
end program print_test

max_stars_repo_name: caguerra/Fortran-MOOC

max_issues_repo_name: caguerra/Fortran-MOOC

max_forks_repo_name: caguerra/Fortran-MOOC


! { dg-do run }
!
  character(*), parameter  :: chrs = '-+.0123456789eEdD'
  character(*), parameter  :: expr = '-+.0123456789eEdD'
  integer :: i

  if (index(chrs(:), expr) /= 1) STOP 1
  if (index(chrs(14:), expr) /= 0) STOP 2
  if (index(chrs(:12), expr) /= 0) STOP 3
  if (index(chrs, expr(:)) /= 1) STOP 4
  if (index(chrs, expr(1:)) /= 1) STOP 5
  if (index(chrs, expr(:1)) /= 1) STOP 6

  if (foo(expr) /= 1) STOP 7
  if (foo(expr) /= 1) STOP 8
  if (foo(expr) /= 1) STOP 9
  if (foo(expr(:)) /= 1) STOP 10
  if (foo(expr(1:)) /= 1) STOP 11
  if (foo(expr(:1)) /= 1) STOP 12

  call bar(expr)

contains
  subroutine bar(expr)
    character(*), intent(in) :: expr
    character(*), parameter  :: chrs = '-+.0123456789eEdD'
    integer :: foo

    if (index(chrs(:), expr) /= 1) STOP 13
    if (index(chrs(14:), expr) /= 0) STOP 14
    if (index(chrs(:12), expr) /= 0) STOP 15
    if (index(chrs, expr(:)) /= 1) STOP 16
    if (index(chrs, expr(1:)) /= 1) STOP 17
    if (index(chrs, expr(:1)) /= 1) STOP 18
  end subroutine bar

  integer function foo(expr)
    character(*), intent(in) :: expr
    character(*), parameter  :: chrs = '-+.0123456789eEdD'

    foo = index(chrs, expr)
  end function foo

end


text:
!    SpiralMaker software (General Node Model)
!    Computational model to simulate embryonic cleavage.
!    Copyright (C) 2014 Miquel Marin-Riera, Miguel Brun-Usan, Roland Zimm, Tommi Välikangas & Isaac Salazar-Ciudad

!    This program is free software: you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation, either version 3 of the License, or
!    any later version.

!    This program is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.

!    You should have received a copy of the GNU General Public License
!    along with this program.  If not, see <http://www.gnu.org/licenses/>.

module biomechanic	!>>>>>>>> by Miquel 17-5-13

use general
use genetic
use neighboring

contains

subroutine iterdiferencial
integer::nodmo,i,j,k,ii
real*8::a,b,c
real*8::ox,oy,oz !miguel4-1-13

  !CALCULATING FORCES AND MOVEMENT VECTORS

  call forces

  a=0.0d0
  do i=1,nd
    if (dex(i)>a) then ; a=dex(i) ; ii=i ; end if
  end do

  if (a<epsilod) then ; rtime=rtime+delta ; delta=deltamin; return ; end if  !>>> Is 29-8-13
  if(ffu(12)==0)then
    delta=resmax/a      !delta is adjusted so the node that has the greatest force
                      !will always move the same absolute distance, and the rest will
                      !move proportionally to that force
  else
    delta=deltamin
!    delta=deltamax
  end if

  if (delta>deltamax) delta=deltamax
  if (delta<deltamin) delta=deltamin

end subroutine iterdiferencial

!***************************************************************************************************

subroutine rungekutta4(d)  ! Runge-Kutta fourth order integration
real*8 d,halfd,sixthd
real*8 ox(nd),oy(nd),oz(nd)
real*8 kux(nd),kuy(nd),kuz(nd)
real*8 kdx(nd),kdy(nd),kdz(nd)
real*8 ktx(nd),kty(nd),ktz(nd)
real*8 kqx(nd),kqy(nd),kqz(nd)

halfd=d*0.5d0
sixthd=d/6.0d0

ox=node(:nd)%x ; oy=node(:nd)%y ; oz=node(:nd)%z

!k1
kux=px(:nd) ; kuy=py(:nd) ; kuz=pz(:nd)

!k2
node(:nd)%x=node(:nd)%x+halfd*px(:nd)
node(:nd)%y=node(:nd)%y+halfd*py(:nd)
node(:nd)%z=node(:nd)%z+halfd*pz(:nd)

if(nd>1) call neighbor_build
call forces

kdx=px(:nd) ; kdy=py(:nd) ; kdz=pz(:nd)

!k3
node(:nd)%x=node(:nd)%x+halfd*px(:nd)
node(:nd)%y=node(:nd)%y+halfd*py(:nd)
node(:nd)%z=node(:nd)%z+halfd*pz(:nd)

if(nd>1) call neighbor_build
call forces

ktx=px(:nd) ; kty=py(:nd) ; ktz=pz(:nd)

!k4
node(:nd)%x=node(:nd)%x+d*px(:nd)
node(:nd)%y=node(:nd)%y+d*py(:nd)
node(:nd)%z=node(:nd)%z+d*pz(:nd)

if(nd>1) call neighbor_build
call forces

kqx=px(:nd) ; kqy=py(:nd) ; kqz=pz(:nd)

!final
node(:nd)%x=ox+sixthd*(kux+2*kdx+2*ktx+kqx)
node(:nd)%y=oy+sixthd*(kuy+2*kdy+2*kty+kqy)
node(:nd)%z=oz+sixthd*(kuz+2*kdz+2*ktz+kqz)

end subroutine

!***************************************************************************************************

subroutine adaptive_rungekutta
real*8 ox(nd),oy(nd),oz(nd)
real*8 aux(nd),auy(nd),auz(nd)
real*8 adx(nd),ady(nd),adz(nd)
real*8 r,halfdelta,invdelta,mdx,mdy,mdz,suggesteddelta

37 continue

halfdelta=0.5d0*delta
invdelta=1d0/delta

ox=node(:nd)%x ; oy=node(:nd)%y ; oz=node(:nd)%z

call rungekutta4(delta)

aux=node(:nd)%x ; auy=node(:nd)%y ; auz=node(:nd)%z
node(:nd)%x=ox  ; node(:nd)%y=oy  ; node(:nd)%z=oz

call rungekutta4(halfdelta)
call rungekutta4(halfdelta)

adx=node(:nd)%x ; ady=node(:nd)%y ; adz=node(:nd)%z

mdx=maxval(abs(aux-adx)*invdelta)
mdy=maxval(abs(auy-ady)*invdelta)
mdz=maxval(abs(auz-adz)*invdelta)

if (mdx>=mdy.and.mdx>=mdz) r=mdx
if (mdy>=mdz.and.mdy>=mdz) r=mdy
if (mdz>=mdy.and.mdz>=mdx) r=mdz

suggesteddelta=0.9d0*prec*delta/r

if (r>prec) then !the step is too large
  delta=suggesteddelta
!print *,"NO",delta,r,prec
  goto 37
else
!print *,"SI",delta,r,prec
  ! in theory we should make delta=suggested delta but we prefer delta to be decided based on resmax in each step
  node(:nd)%x=adx
  node(:nd)%y=ady
  node(:nd)%z=adz
end if

end subroutine

!***************************************************************************************************

subroutine forces
real*8   ::ix,iy,iz,dd
real*8   ::a,b,c,d,e,f,g
integer  ::celi,celj,nod
real*8   ::youe,repe,adhe,adho,repcele,deqe,ideqe
real*8   ::younod,repnod,adhnod,repcelnod,reqnod,tornod,stornod    !>>>> Miquel 16-8-13
real*8   ::ax,ay,az,bx,by,bz
real*8   ::ud,udd,uddd
real*8   ::cx,cy,cz,ccx,ccy,ccz,dotp,pesco
real*8   ::icx,icy,icz,idd,iudd,id !>>>>>>>> MIQUEL 4-3-13
real*8   ::mcx,mcy,mcz,md,umd					!>>>>>>>> MIQUEL 30-4-13
real*8   ::nodda,posca
integer  ::ivv		!>>>>>>>> MIQUEL 4-3-13
integer  ::i,j,ii,jj,kk,ic,iii,jjj,kkk,iiii,jjjj,kkkk,iv,kjjj,jkkk
integer  ::nuve,inuve
integer  ::tipi,tipic																!>>>>>>>>>>>>>>>>>>>>>Miquel 23-4-13
integer  ::switch,twoep

integer  ::lateral,vertical !flags that tell if there is a lateral or vertical component to take into account !>>Miquel28-1-14

real*8   ::rvx,rvy,rvz   !the resulting force vector
real*8   ::uvx,uvy,uvz   !unit vector
real*8   ::pox,poy,poz   !polarisation vector (from the cell)

real*8   ::ad,fd  !>>Miquel28-1-14

real*8   ::upr(nd)
integer  ::iupr(nd)
real*8   ::r(nd),er(nd)

real*8   ::rcilx(nd),rcily(nd),rcilz(nd)
real*8   ::rtorx(nd),rtory(nd),rtorz(nd)
real*8   ::rstorx(nd),rstory(nd),rstorz(nd)
real*8   ::rsprx,rspry,rsprz

!integer, parameter :: mnn=200  !declared in general module !>>Miquel24-2-14

                              ! ACHTUNG, WE ASSUME THAT THE MAXIMAL NUMBER OF NODES INTERACTING WITH A NODE IS mnn, otherwise these
                              ! matrices become unbareably large for large systems
!integer  ::vei(omnn,nd)        ! force storing arrays for optimization    >>>>>>Miquel 7-8-16
real*8   ::vforce(omnn,nd)     ! so we don't calculate each interaction twice
real*8   ::vuvx(omnn,nd)       ! (for now we do it only for cylinder and normal interactions,
real*8   ::vuvy(omnn,nd)       ! not torsion nor springs)
real*8   ::vuvz(omnn,nd)       !
real*8   ::vtorforce(omnn,nd)  ! NOTE THE NON-CONVENTIONAL FORM OF THE ARRAYS:
real*8   ::vtoruvx(omnn,nd)    ! DIM 1 IS THE LISTS THE NEIGHBORS,
real*8   ::vtoruvy(omnn,nd)    ! DIM 2 LISTS THE NODE ID
real*8   ::vtoruvz(omnn,nd)    !
real*8   ::vstorforce(omnn,nd) !
real*8   ::vstoruvx(omnn,nd)   !
real*8   ::vstoruvy(omnn,nd)   !
real*8   ::vstoruvz(omnn,nd)   !
!integer ::nveins(nd)         !
integer  ::epinveins(nd)      ! we store how many same-side epithelial neighbors an epithelial node has !>>>Miquel4-4-14
integer  ::alone              ! 0 if the node is really alone

  vcilx=0 ; vcily=0 ; vcilz=0 ; vsprx=0     !force vectors storage matrices, for different components
  vtorx=0 ; vtory=0 ; vtorz=0 ; vspry=0     !of the resulting force
  vstorx=0 ; vstory=0 ; vstorz=0 ; vsprz=0
  vei=0 ; vforce=0 ; vuvx=0 ; vuvy=0 ; vuvz=0 ; nveins=0  !>>>>>>Miquel 7-8-16
  vtorforce=0 ; vtoruvx=0 ; vtoruvy=0 ; vtoruvz=0         !
  vstorforce=0 ; vstoruvx=0 ; vstoruvy=0 ; vstoruvz=0     !

  fmeanl=0 ; fmeanv=0  !storage vector that makes the balance between compressive and tensile forces within a node    !>>>Miquel23-1-14

  rcilx=0d0  ; rcily=0d0  ; rcilz=0d0  !they store the force components for all the nodes !>>Miquel4-4-14
  rtorx=0d0  ; rtory=0d0  ; rtorz=0d0
  rstorx=0d0 ; rstory=0d0 ; rstorz=0d0
  epinveins=0

  do nod=1,nd
    !lonely=0 !fossile? >>Miquel27-2-14

    if (node(nod)%hold==2) then                ! >>> Is 30-6-14
      dex(nod)=0  !module of the force vector ! >>> Is 30-6-14
      px(nod)=0 ; py(nod)=0 ; pz(nod)=0       ! >>> Is 30-6-14
      cycle                                   ! >>> Is 30-6-14
    end if                                    ! >>> Is 30-6-14

    rsprx=0.0d0  ; rspry=0.0d0; rsprz=0.0d0
    ix=node(nod)%x ; iy=node(nod)%y ; iz=node(nod)%z
    tipi=node(nod)%tipus ; celi=node(nod)%icel
    iii=nint(ix*urv)    ; jjj=nint(iy*urv)   ; kkk=nint(iz*urv)
    rvx=0d0    ; rvy=0d0    ; rvz=0d0
    nuve=0!nuve=nveins(nod) !>>>>>Miquel 7-8-13 : this is not always zero since we fill this matrix from its neighbours
    switch=0
    alone=0

    !SPRINGS

      iv=0    !>>>>Miquel17-1-14

    younod=node(nod)%you                                            !>>>>>Miquel 16-8-13
    repnod=node(nod)%rep                                            !
    adhnod=node(nod)%adh !default inespecific adhesion of node nodmo!
    repcelnod=node(nod)%repcel                                      !
    tornod=node(nod)%tor                                            !
    stornod=node(nod)%stor                                          !
    reqnod=node(nod)%req                                            !
    nodda=node(nod)%da

    !NODE's REPULSIONS AND ADHESIONS

    do i=1,nneigh(nod)
      ic=neigh(nod,i)
      if(ic<nod)then !we have calculated that interaction already !>>Miquel4-4-14
        alone=1
        cycle
      end if

      !so it turns out that the nod-ic interactions has not been calculated before

      bx=node(ic)%x   ; by=node(ic)%y    ; bz=node(ic)%z
      ccx=bx-ix       ; ccy=by-iy        ; ccz=bz-iz
      !d=sqrt(ccx**2+ccy**2+ccz**2)
      d=dneigh(nod,i)
      ud=1d0/d
      tipic=node(ic)%tipus
      twoep=0
      ad=0 ; fd=0 ; ddd=0     !>>Miquel28-1-14

      alone=1      !this is crappy but fast, it makes that lonely nodes are eliminated in squares

          fd=d !BOTH NODES ARE NON-EPITHELIAL: we just consider the interactions between nodes as such
          if (fd-nodda-node(ic)%da>epsilod) cycle
          uvx=ccx*ud ; uvy=ccy*ud ; uvz=ccz*ud  !unit vector of the within cilinder force !el modul és el mateix que el vector c

300   nuve=nuve+1

      if (nuve>mnn) then;
        print *,"PANIC!!! PANIC!!!PANIC!!!; this is going to crash because there is too many " ;
        print *,"neighbors per a node to interact, and vuvx matrix and those are too small"
      end if

      !ALL THAT WAS JUST TO CALCULATE THE RIGHT DISTANCE BETWEEN NODES, ddd, NOW WE CALCULATE THE ACTUAL ENERGIES
      !if (ffu(2)==1) then IS 25-12-13
        if(node(nod)%icel==node(ic)%icel)then
          deqe=(reqnod+node(ic)%req)             !>>>> Miquel 16-8-13
          if(fd-deqe<-epsilod)then
            repe=0.5d0*(repnod+node(ic)%rep)     !>>>> Miquel 16-8-13
            f=2*repe*(fd-deqe)
          else
            youe=(younod+node(ic)%you)          !>>>> Miquel 16-8-13
            f=youe*(fd-deqe)
          end if
        else
          deqe=(reqnod+node(ic)%req)          !>>>> Miquel 16-8-13
          if(fd-deqe<-epsilod)then
            repcele=(repcelnod+node(ic)%repcel) !>>>> Miquel 16-8-13
            f=repcele*(fd-deqe)       !in fact that is the derivative of the energy
          else
            adhe=0.5d0*(adhnod+node(ic)%adh)          !>>>> Miquel 16-8-13
            if (npag(1)>0) then ! we have adhesion molecules
              do j=1,npag(1)
                k=whonpag(1,j)
                do kjjj=1,npag(1)
                  jkkk=whonpag(1,kjjj)
                  if (gex(nod,k)>0.0d0.and.gex(ic,jkkk)>0.0d0) then
                    adhe=adhe+gex(nod,k)*gex(ic,jkkk)*kadh(int(gen(k)%wa(1)),int(gen(jkkk)%wa(1))) ! >>> Is 7-6-14    !this is specific adhesion
                  end if
                end do
              end do
            end if
            f=2*adhe*(fd-deqe)          !in fact that is the derivative of the energy
          end if
        end if

        rcilx(nod)=rcilx(nod)+f*uvx ; rcily(nod)=rcily(nod)+f*uvy ; rcilz(nod)=rcilz(nod)+f*uvz  !>>>Miquel 7-8-13
        rcilx(ic)=rcilx(ic)-f*uvx ; rcily(ic)=rcily(ic)-f*uvy ; rcilz(ic)=rcilz(ic)-f*uvz  !>>>Miquel 4-4-14

        fmeanl(nod)=fmeanl(nod)+(fd-deqe) ; fmeanl(ic)=fmeanl(ic)+(fd-deqe) ! >>> Is 21-6-14

324   continue
    end do

    if(epinveins(nod)>0)then                  !>>>Miquel24-3-14
      fmeanl(nod)=fmeanl(nod)/epinveins(nod)
    end if

789 if (ffu(8)==1) then   !this is kind of crappy in the sense that
      if (alone==0) then
        node(nod)%talone=node(nod)%talone+1
      else
        node(nod)%talone=0
      end if
    end if

    if (aut==0) then ! if aut==0 there is not display
      vsprx(nod)=rsprx ; vspry(nod)=rspry ; vsprz(nod)=rsprz                     !putting the force components into storage vectors for display
      vcilx(nod)=rcilx(nod)    ; vcily(nod)=rcily(nod)    ; vcilz(nod)=rcilz(nod)         !this part should be turned down
      vtorx(nod)=rtorx(nod)    ; vtory(nod)=rtory(nod)    ; vtorz(nod)=rtorz(nod)         !when there is no display
      vstorx(nod)=rstorx(nod)  ; vstory(nod)=rstory(nod)  ; vstorz(nod)=rstorz(nod)!
    end if

    if (epinveins(nod)>0) then  ! >>> Is 7-6-14
      a=epinveins(nod) ! >>> Is 7-6-14
      a=1.0d0/a        ! >>> Is 7-6-14
      rvx=rsprx+rcilx(nod)+(rtorx(nod)+rstorx(nod))*a  ! >>> Is 7-6-14 !summing the force components into a net force vector
      rvy=rspry+rcily(nod)+(rtory(nod)+rstory(nod))*a  ! >>> Is 7-6-14
      rvz=rsprz+rcilz(nod)+(rtorz(nod)+rstorz(nod))*a  ! >>> Is 7-6-14
    else
      rvx=rsprx+rcilx(nod)  ! >>> Is 7-6-14 !summing the force components into a net force vector
      rvy=rspry+rcily(nod)  ! >>> Is 7-6-14
      rvz=rsprz+rcilz(nod)  ! >>> Is 7-6-14
    end if

    !Physical boundaries force  !>>>>>Miquel9-1-14
    if(node(nod)%hold==1)then
      uvx=node(nod)%orix-ix
      uvy=node(nod)%oriy-iy
      uvz=node(nod)%oriz-iz
      d=uvx**2+uvy**2+uvz**2 ; if(d>epsilod)then ; ud=1d0/sqrt(d);else;d=0;end if
      a=node(nod)%khold
      rvx=rvx+uvx*a*ud    !now the force is constant        !no need to calculate the unit vector, because here you have the product of the unit vector and the distance,
      rvy=rvy+uvy*a*ud    !to make it a spring remove ud    !wich is the original vector
      rvz=rvz+uvz*a*ud
    end if

    dex(nod)=sqrt(rvx**2+rvy**2+rvz**2)  !module of the force vector
    px(nod)=rvx ; py(nod)=rvy ; pz(nod)=rvz
  end do

end subroutine forces

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine ordenarepeq(ma,mt,rang)
  integer rang
  real*8 ma(rang)
  integer mt(rang)
  integer i,j,k
  real*8 a
    mt=0
el: do i=1,rang
      a=ma(i) ; k=1
      do j=1,rang ; if (a>ma(j)) k=k+1 ; end do
      do j=k,rang ; if (mt(j)==0) then ; mt(j)=i ; cycle el ; end if ; end do
    end do el
end subroutine ordenarepeq

end module biomechanic

max_stars_repo_name: isalazarciudad/SpiralMaker

max_issues_repo_name: sebjameswml/SpiralMaker

max_forks_repo_name: sebjameswml/SpiralMaker


text:
      subroutine adecay(p,pe,pnb,pc,m)
      implicit none
************************************************************************
*     Author: R.K. Ellis, January 2012                                 *
*     antitop decay  a --> e(pe)+nb(pnb)+bbar(pc)                      *
*     with bottom and top masses (and no radiation)                    *
*     in massless spinor notation                                      *
*     pe,pnb,pc are integers that point to                             *
*     the appropriate four-momenta in p                                *
*     pe=electron                                                      *
*     pnb=antineutrino                                                 *
*     pc=anti-bottom quark                                             *
*     q(c) is rendered massless wrt to pnb                             *
*     q(a) is rendered massless wrt to pe                              *
*     returned m(apol,cpol)                                            *
************************************************************************
      include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      include 'masses.f'
      double precision p(mxpart,4),q(mxpart,4),dot,sw,ala,alc
      double complex m(2,2),cprop
      integer nb,a,c,e,si,pc,pe,pnb
      parameter(a=1,e=3,nb=4,c=2)
      do si=1,4
      q(a,si)=p(pe,si)+p(pnb,si)+p(pc,si)
      q(e,si)=p(pe,si)
      q(nb,si)=p(pnb,si)
      q(c,si)=p(pc,si)
      enddo
      ala=mt**2/(2d0*dot(q,a,e))
      alc=mb**2/(2d0*dot(q,c,nb))
      do si=1,4
      q(a,si)=q(a,si)-ala*q(e,si)
      q(c,si)=q(c,si)-alc*q(nb,si)
      enddo

      call spinoru(4,q,za,zb)
      sw=s(e,nb)
      cprop=dcmplx(sw-wmass**2,wmass*wwidth)
C---order of polarizations is the m(apol,cpol)
      m(1,1)=czip
      m(1,2)=czip
      m(2,1)=czip
      m(2,2)=-za(e,a)*zb(c,nb)/cprop
      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
SUBROUTINE gthr

!-------------------------------------------------------------
!
! Gathers all the f vectors and puts solution into flux
!
!-------------------------------------------------------------

USE invar
USE solvar
IMPLICIT NONE
INTEGER :: tag, indx, pi, eqs, i, g, ierr
INTEGER, DIMENSION(3) ::  istat
REAL*8, DIMENSION(:,:), ALLOCATABLE :: tmpf

INCLUDE 'mpif.h'

! Bring all the blocks' solutions for all the groups back to root for output
IF (irank /= root) THEN
   tag = 100 + irank
   CALL MPI_SEND(f,nx*sord*ng,MPI_DOUBLE_PRECISION,root,tag,MPI_COMM_WORLD,ierr)
ELSE
   indx = 0
   DO pi = 0, npx-1
      eqs = nxvec(pi)*sord*ng
      IF (pi /= 0) THEN
         ALLOCATE(tmpf(nxvec(pi)*sord,ng))
         tmpf = 0.0
         tag = 100 + pi
         CALL MPI_RECV(tmpf,eqs,MPI_DOUBLE_PRECISION,pi,tag,MPI_COMM_WORLD,istat,ierr)
         DO g = 1, ng
            DO i = 1, nxvec(pi)*sord
               flux(i+indx,g) = tmpf(i,g)
            END DO
         END DO
         DEALLOCATE(tmpf)
      ELSE IF (pi == 0) THEN
         ! Root copies its own solution into flux
         DO g = 1, ng
            DO i = 1, nx*sord
               flux(i+indx,g) = f(i,g)
            END DO
         END DO
      END IF
   indx = indx + nxvec(pi)*sord
   END DO
END IF

IF (irank == root) DEALLOCATE(nxvec)

RETURN
END SUBROUTINE gthr

max_stars_repo_name: lanl/pidots

max_issues_repo_name: lanl/pidots

max_forks_repo_name: lanl/pidots


text:
!
! Fortran string class
!
! Allows for assignment from character variables and the ability to get /
! set the internal character value.
!
module string_m
    ! Set defaults for class
    implicit none
    private

    ! Declare the type
    type, public :: string_t
        private

        ! Internal private character store
        character(:), allocatable :: str_m
    contains
        private

        ! Getter / setter functions
        procedure, public, pass(this) :: get_value => get_value_string_t
        procedure, public, pass(this) :: set_value => set_value_string_t

        ! Assignment operator for other strings and characters
        generic, public  :: assignment (=) => string_t_assign_string_t, &
                                              string_t_assign_character

        ! Procedure declaration for internal methods
        procedure, private, pass(lhs) :: string_t_assign_string_t, &
                                         string_t_assign_character
    end type string_t

    ! Class constructor
    interface string_t
        module procedure string_t_constructor
    end interface string_t

contains
    !
    ! Constructor - initialize internal storage
    !
    type (string_t) function string_t_constructor()
        string_t_constructor%str_m = ""
    end function string_t_constructor

    !
    ! Assignment - for other variables of type(string_t)
    !
    subroutine string_t_assign_string_t(lhs, rhs)
        class (string_t), intent (inout) :: lhs
        class (string_t), intent (in) :: rhs

        lhs%str_m = rhs%str_m
    end subroutine string_t_assign_string_t

    !
    ! Assignment - for characters and character arrays
    !
    elemental subroutine string_t_assign_character(lhs, rhs)
        class (string_t), intent (inout) :: lhs
        character(len=*), intent (in) :: rhs

        lhs%str_m = rhs
    end subroutine string_t_assign_character

    ! Getter function
    pure function get_value_string_t(this) result(res)
        character (:), allocatable :: res
        class (string_t), intent (in) :: this

        res = this%str_m
    end function get_value_string_t

    ! Setter subroutine
    pure subroutine set_value_string_t(this, the_value)
        class (string_t), intent (inout) :: this
        character(len=*), intent (in) :: the_value

        this%str_m = trim(the_value)
    end subroutine set_value_string_t

end module string_m

max_stars_repo_name: bceverly/fortranString

max_issues_repo_name: bceverly/fortranString

max_forks_repo_name: bceverly/fortranString


text:
      SUBROUTINE MATMPY(A,B,C,N,M)
C
C                  'MATMPY' MULTIPLIES TWO SQUARE MATRICES; A*B = C
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      DIMENSION A(N,N),B(N,N),C(N,N)
C
      DO 10 I=1,N
      DO 10 J=1,N
      C(I,J)=0.D0
      DO 10 K=1,N
   10 C(I,J)=C(I,J) + A(I,K)*B(K,J)
      RETURN
      END

ext: for


text:
c Identical to read7, except for the fact that the cross sections can
c be multiplied by factors entered in input
      implicit real*4 (a-h,o-z)
      real * 4 xfacts(1:4)
      logical flag,titfnd
      character * 70 file,file1,filename(1:4)
      character * 70 line,vartype,vartype2
      character * 70 strscale,strlimx,strlimy,strtit
      character * 11 chref1,chref7
      character * 14 chref2,chref3
      character * 18 chref4,chref6
      character * 7 chref5,chref8
      character * 9 chref9,chref10
      character * 10 chref11
c
999   format(3(d16.6))
      chref1='  SET SCALE '
      chref2='  SET LIMITS X'
      chref3='  SET LIMITS Y'
      chref4='  SET ORDER X Y DY'
      chref5='   PLOT'
      chref6='  ( SET FONT DUPLE'
      chref7='   NEW PLOT'
      chref8='  HIST '
      chref9=' new plot'
      chref10=' (       '
      chref11='  TITLE   '
c
      strtit(1:17)='  TITLE    TOP   '
      write(6,*)'enter number of files (less than 5)'
      read(5,*)nfil
      if(nfil.ge.5)stop
      write(6,*)'enter filenames'
      read(5,'(a)')file
      filename(1)=file
      open(unit=50,type='old',name=file,form='formatted')
      do i=1,nfil-1
        iunit=50+i
        read(5,'(a)')file1
        filename(i+1)=file1
        open(unit=iunit,type='old',name=file1,form='formatted')
      enddo
      write(6,*)'enter multiplicative factors'
      do i=1,nfil
        read(5,*)xfacts(i)
      enddo
c
      iend=0
      imtfill=0
      vartype2='&&&&&&&&'
 111  continue
      titfnd=.false.
      write(88,'(a)')' ('
      write(88,'(a)')' new plot'
      write(88,'(a)')' ('
      flag=.true.
      do while(flag)
        if(imtfill.eq.1)then
          flag=.false.
          write(88,'(a)')strtit
          write(88,'(a)')strscale
          write(88,'(a)')strlimx
          write(88,'(a)')strlimy
          vartype=vartype2
          write(88,'(a)')vartype
          read(50,'(a)',err=200,end=200)line
          if(line(1:7).ne.' ( INT=')then
            write(6,*)'Bad format in input file'
            stop
          else
            write(88,'(a)')line
          endif
        else
          read(50,'(a)',err=200,end=200)line
          if(line(1:10) .eq. chref11 .and. .not.titfnd) then
             istrt=index(line,'"')
             istp =index(line(istrt+1:70),'"')
             if (istrt.ne.0 .and. istp.ne.0) then
                write(strtit(18:18+istp),'(a)')
     &               line(istrt:istp+istrt)
                do i=19+istp,70
                   write(strtit(i:i),'(a)')' '
                enddo
                write(88,'(a)')strtit
                titfnd=.true.
             endif
          endif
          if(line(1:11).eq.chref1)then
            write(88,'(a)')line
            strscale=line
          elseif(line(1:14).eq.chref2)then
            write(88,'(a)')line
            strlimx=line
          elseif(line(1:14).eq.chref3)then
            write(88,'(a)')line
            strlimy=line
          elseif(line(1:18).eq.chref4)then
            flag=.false.
            write(88,'(a)')line
            read(50,'(a)',err=200,end=200)line
            vartype=line
            write(88,'(a)')line
            read(50,'(a)',err=200,end=200)line
            if(line(1:7).ne.' ( INT=')then
              write(6,*)'Bad format in input file'
              stop
            else
              write(88,'(a)')line
            endif
          else
            continue
          endif
        endif
      enddo
      write(99,*)"Doing ",vartype
c now read the numbers
 222  continue
      do jj=1,100000
        read(50,'(a)',err=200,end=200)line
        read(unit=line,fmt=*,err=10)x,y,dy
        write(88,999)x,y*xfacts(1),dy
      enddo
 10   continue
      if(line(1:7).ne.chref8)then
        write(6,*)'No HIST statement at the end of the histogram'
        stop
      else
        write(88,'(a)')'  hist solid'
      endif
      read(50,'(a)',err=200,end=200)line
      if(line(1:7).ne.chref5)then
        write(6,*)'No PLOT statement at the end of the histogram'
        stop
      endif
      read(50,'(a)',err=200,end=333)line
      if(line(1:18).eq.chref6.or.line(1:11).eq.chref7.or.
     #   line(1:9).eq.chref9.or.line(1:9).eq.chref10)then
        imtfill=0
        vartype2='&&&&&&&&'
      else
        imtfill=1
        vartype2=line
      endif
      goto 444
 333  continue
      iend=1
 444  continue
c now the other files
      do i=1,nfil-1
        iunit=50+i
        flag=.true.
        do while(flag)
          read(iunit,'(a)',err=200,end=200)line
          if(line.eq.vartype)then
            flag=.false.
            write(88,'(a)')line
            read(iunit,'(a)',err=200,end=200)line
            if(line(1:7).ne.' ( INT=')then
              write(6,*)'Bad format in input file'
              stop
            else
              write(88,'(a)')line
            endif
          else
            continue
          endif
        enddo
c now read the numbers
        do jj=1,100000
          read(iunit,'(a)',err=200,end=200)line
          read(unit=line,fmt=*,err=11)x,y,dy
          write(88,999)x,y*xfacts(i+1),dy
        enddo
 11     continue
        if(line(1:7).ne.chref8)then
          write(6,*)'No HIST statement at the end of the histogram'
          stop
        else
          if(i.eq.1)then
            write(88,'(a)')' set pattern .05 .07'
          elseif(i.eq.2)then
            write(88,'(a)')' set pattern .01 .07'
          elseif(i.eq.3)then
            write(88,'(a)')' set pattern .01 .07 .05 .07'
          else
            write(6,*)'Fatal error'
            stop
          endif
          write(88,'(a)')' hist patterned'
        endif
        if(iend.eq.0)then
          read(iunit,'(a)',err=200,end=200)line
          if(line(1:7).ne.chref5)then
            write(6,*)'No PLOT statement at the end of the histogram'
            stop
          endif
        endif
      enddo
c
      if(iend.eq.0)goto 111
200   continue
      end

      function istrl(string)
c returns the position of the last non-blank character in string
      character * (*) string
      i = len(string)
      dowhile(i.gt.0.and.string(i:i).eq.' ')
         i=i-1
      enddo
      istrl = i
      end

      subroutine strcat(str1,str2,str)
c concatenates str1 and str2 into str. Ignores trailing blanks of str1,str2
      character *(*) str1,str2,str
      l1=istrl(str1)
      l2=istrl(str2)
      l =len(str)
      if(l.lt.l1+l2) then
          write(*,*) 'error: l1+l2>l in strcat'
          write(*,*) 'l1=',l1,' str1=',str1
          write(*,*) 'l2=',l2,' str2=',str2
          write(*,*) 'l=',l
          stop
      endif
      if(l1.ne.0) str(1:l1)=str1(1:l1)
      if(l2.ne.0) str(l1+1:l1+l2)=str2(1:l2)
      if(l1+l2+1.le.l) str(l1+l2+1:l)= ' '
      end

max_stars_repo_name: JaySandesara/madminer

max_issues_repo_name: JaySandesara/madminer

max_forks_repo_name: JaySandesara/madminer


text:
PROGRAM read_n_calc
IMPLICIT NONE
REAL, ALLOCATABLE, DIMENSION(:,:) :: data
CHARACTER(LEN=20) :: filename
INTEGER :: status
INTEGER :: rows, cols, i, j
CHARACTER(LEN=10) :: row_str, col_str
CHARACTER(LEN=40) :: row_sum_fmt, col_sum_fmt

WRITE(*, *) "Enter the filename to read and process:"
READ(*, *) filename

OPEN(1, FILE=filename, STATUS='OLD', ACTION='READ', IOSTAT=status)
READ(1, *) rows, cols
ALLOCATE(data(1:rows, 1:cols))

WRITE(row_str, '(I3)') rows
WRITE(col_str, '(I3)') cols
row_str = TRIM(row_str)
col_str = TRIM(col_str)
row_sum_fmt = '(' // row_str // '(A, I4, 2X, A, F12.1, A))'
col_sum_fmt = '(' // col_str // '(A, I4, 2X, A, F12.1, A))'
row_sum_fmt = TRIM(row_sum_fmt)
col_sum_fmt = TRIM(col_sum_fmt)

DO i = 1, rows
  READ(1, '(' // col_str // 'F8.1)', IOSTAT=status) (data(i, j), j = 1, cols)
  IF (status .lt. 0) THEN
    EXIT
  END IF
END DO

WRITE(*, row_sum_fmt) ("Row:", i, "Sum:", SUM(data(i, 1:cols)), NEW_LINE('A'), i = 1, rows)
WRITE(*, col_sum_fmt) ("Col:", i, "Sum:", SUM(data(1:rows, i)), NEW_LINE('A'), i = 1, cols)
!WRITE(*, row_sum_fmt) ("Row:", i, "Sum:", SUM(data(i, 1:cols)), NEW_LINE('A'), i = 1, rows)
!WRITE(*, col_sum_fmt) ("Col:", i, "Sum:", SUM(data(1:rows, i)), i = 1, cols)

DEALLOCATE(data)
END PROGRAM read_n_calc

max_stars_repo_name: evanmacbride/fortran-practice

max_issues_repo_name: evanmacbride/fortran-practice

max_forks_repo_name: evanmacbride/fortran-practice


text:
      subroutine sol_print(n, sn, x, blo, bup, f, vr_reg, k_it,
     &     lws, nlws, ws, nws)
c
c     sol_print: prints the solution in sol####.dat and sol####.m
c                in subdirectory <dirn>
c
c     f           optimal value
c     n           #vars
c     x(n)        value of vars
c     blo(n)      bounds on vars
c     bup(n)
c     vr_reg      var registering information
c     lws(nlws)   workspace
c     dirn        directory to write file to
c     k_it        #of solve (this is used for the #### in filename)
c
      use interface_mod

      implicit none

      INCLUDE 'SLPCOMM.INC'
      include 'pusr.inc'

      integer n, k_it, nlws, sn, nws
      double precision f

      double precision x(n), blo(n), bup(n), ws(nws)
      integer lws(nlws), vr_reg(2+2*sn+3*n)

      integer i, j, k, l, demin, demax

      character*4 int2char
      integer trim77

      integer        seed, bigiter
      common /summc/ seed, bigiter

        DIMENSION ITEMP(2)
        DOUBLE PRECISION DTEMP
        INTEGER ITEMP
        EQUIVALENCE (ITEMP,DTEMP)

c     ==================== subroutine body =======================

      if (max((mi+de)*(rm+mi), (nu+1)*(mi+de)).gt.nlws) then
         write(nout, *) 'Not enough integer workspace in solprint'
         stop
      end if
      if (mi+de.gt.nws) then
         write(nout, *) 'Not enough double workspace in solprint'
         stop
      end if

      open(12, file='sol'//int2char(k_it)//'_.dat')
      do i=1,n
         write(12,*) i, x(i)
      end do
      close(12)

      open(12,file='sol'//int2char(k_it)//'.dat')

      open(13,file='sol'//int2char(k_it)//'.m')

      write(12,*)  ' ============= SOLUTION ================='
      write(12,*)
      write(12,*)  ' f = ',f
      write(13,*) 'f = ',f
      write(12,*)
      write(12,*)  'Composition of mixers (in terms of raw mat & bins):'
      write(12,*)

c     ... get mu's
      do i=1,(mi+de)*(rm+mi)
         lws(i) = 0
      end do

      do i=1,n
         if (vr_reg(2+2*sn+i).eq.TYPE_MU) then
            j = vr_reg(2+2*sn+n+i)
            k = vr_reg(2+2*sn+2*n+i)
            lws((j-1)*(mi+de)+k) = i
         end if
      end do

      do i=1,rm
         do j=1,mi
            if (lws((i-1)*(mi+de)+j).gt.0) then
               ws(j) = x(lws((i-1)*(mi+de)+j))
            else
               ws(j) = 0.d0
            end if
         end do
         write(12,'(A,I3,400F8.3)') 'Rm-',i,(ws(j), j=1,mi)
      end do
      write(12,*) '---------------------------------------------------'
      do i=1,mi
         do j=1,mi
            if (lws((i+rm-1)*(mi+de)+j).gt.0) then
               ws(j) = x(lws((i+rm-1)*(mi+de)+j))
            else
               ws(j) = 0.d0
            end if
         end do
         write(12,'(A,I3,400F8.3)') 'Mi-',i,(ws(j), j=1,mi)
      end do
      write(12,*)
      write(12,*)  'Composition of demand products (in terms of mi+rm):'
      write(12,*)
      k = 1
 10   continue
      demin = (k-1)*9+1
      demax = min(de,demin+8)
      do i=1,mi
         do j=demin,demax
            if (lws((i+rm-1)*(mi+de)+mi+j).gt.0) then
               ws(j) = x(lws((i+rm-1)*(mi+de)+mi+j))
            else
               ws(j) = 0.d0
            end if
         end do
         write(12,'(A,i3,10F8.3)')
     .             'Mi-',i,(ws(j),j=demin,demax)
      end do
      write(12,*) '---------------------------------------------------'
      do i=1,rm
         do j=demin,demax
            if (lws((i-1)*(mi+de)+mi+j).gt.0) then
               ws(j) = x(lws((i-1)*(mi+de)+mi+j))
            else
               ws(j) = 0.d0
            end if
         end do
         write(12,'(A,i3,10F8.3)')
     .             'RM-',i,(ws(j),j=demin,demax)
      end do
      write(12,*)
      k = k+1
      if (demax.lt.de) goto 10

      write(13,*) 'mu = ['
      do i=1,rm+mi
         do j=1,mi+de
            if (lws((i-1)*(mi+de)+j).gt.0) then
               ws(j) = x(lws((i-1)*(mi+de)+j))
            else
               ws(j) = 0.d0
            end if
         end do
         write(13,'(400F20.10)') (ws(j), j=1,mi+de)
      end do
      write(13,*) '];'
      write(13,*)

c     ... scan for TYPE_M variables

c     ... create a look up table:
c         lws((j-1)*(nu+1)+k     is idx of M(j, k)
c         lws((j-1)*(nu+1)+nu+1  is idx of C(j)
      do i=1,nu*(mi+de)
         lws(i) = 0
      end do
      do i=1,n
         if (vr_reg(2+2*sn+i).eq.TYPE_M) then
            j = vr_reg(2+2*sn+n+i)
            k = vr_reg(2+2*sn+2*n+i)
            lws((j-1)*(nu+1)+k) = i
         end if
         if (vr_reg(2+2*sn+i).eq.TYPE_C) then
            j = vr_reg(2+2*sn+n+i)
c            print *, nu, j, (j-1)*(nu+1)+nu+1
            if (j.ge.1)
     &           lws((j-1)*(nu+1)+nu+1) = i
         end if
      end do

      write(12,*)
      write(12,*)  'Composition of mixers in terms of nutrients'
      do i=1,nu
         write(12,'(A,I3,400F8.3)') 'Nu-',i,(x(lws((j-1)*(nu+1)+i)),
     &        j=1,mi)
      end do
      write(12,*)
      write(12,*)  'Cost of mixers'
      write(12,'(A,10F8.3)') '      ',(x(lws((j-1)*(nu+1)+nu+1)),
     &     j=1,mi)
      write(12,*)
      write(12,*)  'Composition of demand products by nutrients'
      do k=1,de
         write(12,'(A,I3,A,F10.5)') 'demand product ',k,
     .                       ' Cost =',x(lws((k+mi-1)*(nu+1)+nu+1))
         do l=1,nu
            write(12,'(A,I3,3F10.5)') 'Nu-',l,
     .       blo(lws((k+mi-1)*(nu+1)+l)),x(lws((k+mi-1)*(nu+1)+l)),
     .       bup(lws((k+mi-1)*(nu+1)+l))
         end do
         write(12,*)
      end do

      write(13,*) 'm = ['
      do i=1,mi+de
         write(13,'(400F20.10)') (x(lws((i-1)*(nu+1)+j)), j=1,nu)
      end do
      write(13,*) '];'
      write(13,*)
      write(13,*) 'mlb = ['
      do i=1,mi+de
         write(13,'(400F20.10)') (blo(lws((i-1)*(nu+1)+j)), j=1,nu)
      end do
      write(13,*) '];'
      write(13,*)
      write(13,*) 'mub = ['
      do i=1,mi+de
         write(13,'(400F20.10)') (bup(lws((i-1)*(nu+1)+j)), j=1,nu)
      end do
      write(13,*) '];'
      write(13,*)
      write(13,*) 'c = ['
      write(13,'(400F20.10)') (x(lws((j-1)*(nu+1)+nu+1)), j=1,mi+de)
      write(13,*) '];'
      write(13,*)

      WRITE(nout,*) ' '
      WRITE(nout,*) ' write solution to file: sol'//int2char(k_it)//'.m'
      WRITE(nout,*) ' '

CJCH      open(12,file=dirname(1:trim77(dirname))//'/sol'//
CJCH     .                                  int2char(k_it)//'.m')

c      do j=1,mi
c         ws(j) = 0.
c         do k=1,de
c            ws(j) = ws(j) + x(pv_lamm+(k-1)*mi+j)*user(rm+nu*rm+k)
c         end do
c      end do
c      write(12,*) 'tonmi= ['
c      do i=1,mi
c         write(12,'(F20.10)') ws(i)
c      end do
c      write(12,*) '];'
c      write(12,*)
c      write(12,*) 'totnumi= ['
c      do i=1,nu
c         write(12,'(20F20.10)') (x(pv_m+(j-1)*nu+i)*ws(j), j=1,mi)
c      end do
c      write(12,*) '];'
c      close(12)

C----------------------- CODE INSERTED BY JCH ------------------------
c
c100     DO 140 J=1,DE
c        WRITE(20,'(I5,A,A)') J,COMMA,'Product'
c           DO 120 I=1,RM
c           DTEMP=x(pv_lamr+(j-1)*rm+i)
c           WRITE(20,'(2Z8,1X,G16.8,A,I5,A,A)')
c     &                                  ITEMP,
c     &                                  x(pv_lamr+(j-1)*rm+i),COMMA,
c     &                                  I,COMMA,'Rawmat'
c120        CONTINUE
c           DO 130 I=1,MI
c           DTEMP=x(pv_lamm+(j-1)*mi+i)
c           WRITE(20,'(2Z8,1X,G16.8,A,I5,A,A)')
c     &                                  ITEMP,
c     &                                  x(pv_lamm+(j-1)*mi+i),COMMA,
c     &                                  I,COMMA,'Premix'
c130        CONTINUE
c140     CONTINUE
c
c        DO 160 J=1,MI
c        WRITE(20,'(I5,A,A)') J,COMMA,'Premix'
c           DO 150 I=1,RM
c           DTEMP=x(pv_mu+(j-1)*rm+i)
c           WRITE(20,'(2Z8,1X,G16.8,A,I5,A,A)')
c     &                                  ITEMP,
c     &                                  x(pv_mu+(j-1)*rm+i),COMMA,
c     &                                  I,COMMA,'Rawmat'
c150        CONTINUE
c160     CONTINUE
c
cC----------------------------------------------------------------------

      end

max_stars_repo_name: ERGO-Code/FortranFilterSLP

max_issues_repo_name: ERGO-Code/FortranFilterSLP

max_forks_repo_name: ERGO-Code/FortranFilterSLP


text:
      SUBROUTINE INAP(BUFIN,NSAMPS)
C     INAP RELIEVES THE USER OF ALL THE BULLSHIT IN ASSIGNING AND TRANSFERRING
C   THE INITIAL DATA TO THE AP IN A SUBROUTINE.
C     THE SIO CONVENTION IS TO ALWAYS RELEASE THE AP BACK TO THE SYSTEM AS SOON
C   AS POSSIBLE SO OTHER USERS CAN USE IT, BUT THE PROBLEM WITH THAT IS THAT
C   APINIT HAS AN UNGODLY AMOUNT OF OVERHEAD.
C     IF THERE ISN'T AN AP ON THE SYSTEM, THE DATA IS IN COMMON /APMEM/.
C  THIS MEMORY JUST LOOKS LIKE A REAL AP DATA MEMORY SO THAT ALL THE INDEXES
C  AND DATA ADDRESSES CAN REMAIN THE SAME, REGARDLESS OF AP EXISTENCE.
C
C  ARGUMENTS:
C    BUF    - THE ARRAY IN THE HOST THAT MAY OR MAY NOT BE IN THE AP.
C    NSAMPS - THE NUMBER OF SAMPLES IN THE ARRAY BUF.
C
C   PAUL HENKART, SCRIPPS OCEANOGRAPHY, NOVEMBER 1979
C
c  mdsize is the size of the ap main data memory.  If there isn't really
c  an ap on the system and we are going to use an ap simulator, which
c  uses an internal array specified in contro, make it very large!
c  Change the ap memory size in contro.f - set mdsize to the size of md!
c
      PARAMETER (IFMT=2)
c      COMMON /APMEM/ A(32768)                                           ! contro has this allocated!
      COMMON /APMEM/ A(5000000)                                          ! contro has this allocated!
      COMMON /SIOAP/ IASGND,IRELSE,IN,IOUT,NEXTAD,LAPSIZ,IFREE,IUSEAP,
     *          IDECIM, mdsize
C          IASGND - AP ASSIGNMENT SWITCH
C                 =0,  AP IS NOT ASSIGNED - CALL APINIT
C                 =1,  AP IS ASSIGNED - DO NOT CALL APINIT
C          IRELSE - AP RELEASE SWITCH
C                 =0,  AP SHOULD BE RELEASED BEFORE RETURNING
C                 =1,  AP SHOULD BE LEFT ASSIGNED
C          IN     - A SWITCH INDICATING WHETHER THE DATA IS ALREADY IN THE AP
C                 =0,  THE DATA IS NOT IN THE AP.
C                 =1,  THE DATA IS IN THE AP AT ADDRESS IN
C          IOUT   - A SWITCH INDICATING WHETHER THE DATA SHOULD BE LEFT IN THE AP
C                 =0,  THE DATA SHOULD BE MOVED BACK TO THE HOST
C                 =1,  THE DATA CAN BE LEFT IN THE AP WITHOUT AN APGET
C          NEXTAD - THE NEXT AVAILABLE FREE LOCATION IN THE AP.  NEXTAD IS SET
C                   TO THE NEXT LOCATION AFTER NSAMPS WHENEVER NEXTAD IS LESS
C                   THAN NSAMPS AND NSAMPS OF DATA WAS PUT INTO THE AP.
C          LAPSIZ - THE INTEGER*4 SIZE OF THE MAIN DATA OF THE AP.
C          IFREE  - A SIGNAL INDICATING THAT SOME SUBROUTINE HAS SOMETHING IN
C                   THE AP THAT MUST BE SAVED, THUS PREVENTING THE AP FROM BEING
C                   RELEASED.
C                 =0, NO SUBROUTINE HAS ANYTHING TO SAVE
C                 <>0, SOMETHINGS ARE IN THE AP THAT MUST BE SAVED - DO NOT
C                      RELEASE THE AP.
C          IUSEAP - A SIGNAL INDICATING WHETHER TO USE THE AP OR NOT
C                 =0, DON'T USE THE AP (THERE MIGHT NOT BE ONE!)
C                 =1, USE THE AP
C****
C****
C****
c  mod 18 Nov 93 - Change the buffer allocation when the data is
c     bigger than originally allocated (when the message about overwriting
c     came out.
c  mod 29 May 96 - Change the initial trace size allocation to be twice the
c                  trace length rather than 1.5 * the trace length.
c  mod 8 Nov 06 - Increase initial buffer size to nsamps*3 (see above), so
c                 FFT of reals has some breathing space.
c   mod 19 Dec 06 - Large changes in nsamps still overwrites good stuff
c  1 Apr 09 - Large changes in nsamps still byting me, so allocate 65536 reals
c             (usigned int) to nextad, so in is always 65536 samples.
c
      INTEGER LSIZE
      INTEGER IASGND,IRELSE,IN,IOUT,NEXTAD,LAPSIZ,IFREE,IUSEAP,IDECIM
      INTEGER NSAMPS
      DIMENSION ISIZE(2)
      EQUIVALENCE (ISIZE(1),LSIZE)
      DIMENSION BUFIN(1)
      SAVE nextad1

      IF(IUSEAP.EQ.0) GO TO 1000                                        ! IF NO AP, DO THE SIMULATOR
      IF(IASGND.EQ.1) GO TO 110
      IASGND=1
      CALL APINIT(0,0,ISTAT)                                             !  INITIALIZE THE AP
      IFREE=0                                                            ! IF ANYONE SAVED ANYTHING TOO BAD - APINIT CLEARS MEMORY!
      IN=0                                                              ! THE DATA CAN'T BE IN THE AP!
      NEXTAD=0                                                          ! LIKEWISE FOR AP ALLOCATION
      IF(ISTAT.GT.0) GO TO 100
      PRINT 10,ISTAT
   10 FORMAT(' ***  ERROR  ***  AP ASSIGNMENT FAILED. STATUS=',I4)
      STOP
  100 CONTINUE
      ISIZE(1)=0                                                        ! CLEAR THE UPPER 16 BITS
      CALL SIZEMD(ISIZE(2))                                             ! THE SIZE USES THE SIGN BIT TO GET TO 64K AP WORDS
      LAPSIZ=LSIZE
      IF(LAPSIZ.LE.0) LAPSIZ = mdsize
  110 IF(IDECIM.LE.0) IDECIM=1
c      IF(NEXTAD.EQ.0) NEXTAD=NSAMPS/IDECIM*2                            !  SET THE NEXT AVAILABLE AP ADDRESS
      IF( nextad .EQ. 0 ) nextad = 65536
      IF(IN.NE.0) RETURN                                                  ! THE DATA IS ALREADY IN THE AP
      IF( NSAMPS .GT. NEXTAD ) THEN
          PRINT 120
          nextad = nsamps + nsamps
      ENDIF
  120 FORMAT(' ***  WARNING  ***  THE DATA GOING INTO THE AP ARE',
     *      ' OVERWRITING PART OF THE RESERVED AP MEMORY.')
      IN=1                                                               ! SET THE AP ADDRESS TO LOCATION 1
      JIN=IN                                                            ! CONVERT TO 16 BIT INTEGER
      JSAMPS=NSAMPS                                                     ! CONVERT TO 16 BIT INTEGER
      CALL APPUT(BUFIN,JIN,JSAMPS,IFMT)
      CALL APWD                                                          ! WAIT FOR THE DATA TRANSFER
      RETURN
C****
C****   TAKE CARE OF THE AP SIMULATOR MEMORY
C****
 1000 IF( IASGND .EQ. 0 ) THEN
          IASGND = 1
          DO 1001 I=1, mdsize
 1001     A(I)=0.                                                            ! AND SOME PROCESS MAY COUNT ON THERE BEING ZEROES THERE
         IFREE=0
         IN=0
         NEXTAD=0
         LAPSIZ = mdsize
      ENDIF
      IF( NEXTAD .EQ. 0 ) THEN
c          NEXTAD=NSAMPS+NSAMPS + nsamps
          nextad = 65536
          nextad1 = nextad
      ENDIF
      IF(IN.NE.0) RETURN                                                ! IS THE DATA ALREADY THERE?
      IF( NSAMPS .GT. NEXTAD1 ) THEN
          PRINT 120
          nextad = nsamps + nsamps + nsamps
          nextad1 = nextad
      ENDIF
      DO 1020 I=1,NSAMPS
 1020 A(I)=BUFIN(I)                                                      ! MOVE THE DATA TO THE SIMULATOR
      IN=1

      RETURN
      END

max_stars_repo_name: apthorpe/sioseis

max_issues_repo_name: henkart/sioseis

max_forks_repo_name: henkart/sioseis


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   EVB-QMDFF - RPMD molecular dynamics and rate constant calculations on
!               black-box generated potential energy surfaces
!
!   Copyright (c) 2021 by Julien Steffen (steffen@pctc.uni-kiel.de)
!                         Stefan Grimme (grimme@thch.uni-bonn.de) (QMDFF code)
!
!   Permission is hereby granted, free of charge, to any person obtaining a
!   copy of this software and associated documentation files (the "Software"),
!   to deal in the Software without restriction, including without limitation
!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!   and/or sell copies of the Software, and to permit persons to whom the
!   Software is furnished to do so, subject to the following conditions:
!
!   The above copyright notice and this permission notice shall be included in
!   all copies or substantial portions of the Software.
!
!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
!   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
!   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
!   THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
!   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
!   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
!   DEALINGS IN THE SOFTWARE.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!     subroutine rpmd_read: read in several settings from qmdff.key file
!      for the rpmd.x program
!
!     part of EVB
!
subroutine rpmd_read(rank,dt,dtdump,dt_info,gen_step,equi_step,umbr_step,xi_min,&
          & xi_max,pmf_min,pmf_max,nbins,pmf_minloc,print_poly,ts_file,recross_mpi)
use general
use evb_mod
use debug
implicit none
!     The current MPI rank
integer, intent(out) :: rank
!     the timestep for calculation and write out
real(kind=8), intent(out) ::dt,dtdump
!     timestep for information of user
real(kind=8), intent(out) ::dt_info
!     several step numbers for trajectories
integer, intent(out) ::equi_step,gen_step,umbr_step
!     borders of XI integration area
real(kind=8), intent(out) :: xi_min,xi_max
!     borders of umbrella windows along Xi reaction path
real(kind=8), intent(out) :: pmf_min,pmf_max
!     number of PMF integratiom bins
integer, intent(out) :: nbins
!     method for location of PMF minimum
character(len=10), intent (out) ::pmf_minloc
!     position of printout trajectory
real(kind=8), intent (out) :: print_poly
!     filename of the ts structure (the starting point!)
character(len=80), intent (out) ::ts_file
!     filename of the educts structure (for unimolecular)
character(len=80)::educts_file
character(len=50)::fix_file ! file with numbers of fixed atoms
!     if recrossing shall be parallelized
logical, intent (out)::recross_mpi
!     Loop indices
integer::i,j,k,l,p
!     Needed characters for read in
character(len=20) keyword
character(len=60)::names,names2
character(len=140) record
character(len=140) string
!     local arrays for breaking/forming bonds
character(len=20),dimension(:),allocatable::input_form,input_break
!     for read in of educt atoms
character(len=1)::at_index(200)
!     auxiliary characters for read in of bonds
character(len=1)::act_char
character(len=20)::act_bond
character(len=20)::act_number
logical::active,exist
integer::next
!     status of actual read in
integer::readstat,istat

!     Initialize several parameters
!
!     number of steps per umbrella structure generation trajectory
gen_step = 0
!     number of steps per umbrella equilibration trajectory
equi_step = 0
!     number of steps per umbrella sampling trajectory
umbr_step = 0
!     number of umbrella trajectories per sampling window
umbr_traj = 0
!     temperature for the whole calculations (Kelvin)
kelvin = 0.0d0
!     timestep for the trajectories (fs)
dt = -1.d0
!     if a trajectory shall be printed out, frequency for printout structures
dtdump = 1.d0
!     number of RPMD beads
nbeads = 0
!     number of equivalent reaction pahts (k(T) multiplier)
npaths = 1
!     if sampling trajectories shall be checked for errors
act_check = .true.
!     if recrossing trajectories shall be checked for errors
recross_check = .true.
!     If the loose check option shall be activated (only development!)
loose_check = .false.
!     The umbrella force constant
k_force = -1.d0
!     The lowest umbrella window (with respect to path progress Xi)
umbr_lo = -5.d0
!     The highest umbrella window (with respect to path progress Xi)
umbr_hi = -5.d0
!     The distance between two umbrella windows
umbr_dist = -1.d0
!     The reaction class with respect to force constant etc.
umbr_type = "none"
!     The asymptotic educts separation for umbrella samplings
r_inf = 0.d0
!     The file in which the startign TS trajectory shall be stored
ts_file = " "
!     The lowest reaction path Xi value for PMF calculation
xi_min = 1.d0
!     The highest reaction path Xi value for PMF calculation
xi_max = -1.d0
!     The number of bins for numerical PMF integration
nbins = 0
!     The method for the calculation of PMF surfaces
pmf_method = "NONE"
!     If the recrossing calculations shall be paralellized
recross_mpi = .false.
!     Number of steps for initial recrossing equilibration trajectory
recr_equi = 0
!     Total number of child trajectories
child_tot = 0
!     timesteps between two child spawning events
child_interv = 0
!     Number of child trajectories per spawning interval
child_point = 0
!     timesteps for a single child trajectory
child_evol = 0
!     Method for location of minimum PMF value on path (asympototic or global min.)
pmf_minloc = "ZERO"
!     If a single umbrella trajectory shall be printed to illustrate RPMD
print_poly = -5.d0
!     The energy tolerance in kJ/mol with respect to the TS energy for errors
energy_tol = 250d0
!     If only the structure generation shall be done for test reasons
gen_test = .false.
!     If the debugging mode shall be set on (write out every single structure)
do_debug = .false.

!
!     Start big loop over all keywords and read them in
!
write(15,*)
write(15,*) "Read in all keywords for the program..."
write(15,*)
do i = 1, nkey
   next = 1
   record = keyline(i)
   call gettext (record,keyword,next)
   call upcase (keyword)
   string = record(next:120)
!
!     reset the read status to successful (check it for each keyword)
!
   readstat= 0
!
!     If a traditional force field shall be calculated, activate global flag!
!     Here also read in a start structure like in dynamic.x to determine the
!     trad. FF coordinates etc.
!

   if (keyword(1:20) .eq. 'GENERATE_STEPS ') then
      read(record,*,iostat=readstat) names,gen_step
   else if (keyword(1:20) .eq. 'EQUILIBR_STEPS ') then
      read(record,*,iostat=readstat) names,equi_step
   else if (keyword(1:20) .eq. 'UMBRELLA_STEPS ') then
      read(record,*,iostat=readstat) names,umbr_step
   else if (keyword(1:20) .eq. 'UMBRELLA_TRAJS ') then
      read(record,*,iostat=readstat) names,umbr_traj
   else if (keyword(1:11) .eq. 'DELTAT ') then
      read(record,*,iostat=readstat) names,dt
   else if (keyword(1:11) .eq. 'TEMP ') then
      read(record,*,iostat=readstat) names,kelvin
   else if (keyword(1:20) .eq. 'THERMOSTAT') then
      read(record,*,iostat=readstat) names,thermo
      call upcase(thermo)
   else if (keyword(1:20) .eq. 'BEAD_NUMBER ') then
      read(record,*,iostat=readstat) names,nbeads
   else if (keyword(1:20) .eq. 'NPATHS ') then
      read(record,*,iostat=readstat) names,npaths
   else if (keyword(1:20) .eq. 'MAX_ERROR ') then
      read(record,*,iostat=readstat) names,err_max
   else if (keyword(1:20) .eq. 'NO_CHECK ') then
      act_check=.false.
   else if (keyword(1:20) .eq. 'RECROSS_NOCHECK ') then
      recross_check=.false.
   else if (keyword(1:20) .eq. 'UMBRELLA_BIAS ') then
      read(record,*,iostat=readstat) names,k_force
   else if (keyword(1:20) .eq. 'UMBRELLA_BONDS ') then
      read(record,*,iostat=readstat) names,umbr_lo,umbr_hi
   else if (keyword(1:20) .eq. 'UMBRELLA_DIST ') then
      read(record,*,iostat=readstat) names,umbr_dist
   else if (keyword(1:15) .eq. 'R_INF ') then
      read(record,*,iostat=readstat) names,r_inf
   else if (keyword(1:20) .eq. 'UMBRELLA_TYPE ') then
      read(record,*,iostat=readstat) names,umbr_type
      call upcase(umbr_type)
   else if (keyword(1:20) .eq. 'PMF_XI_RANGE ') then
      read(record,*,iostat=readstat) names,xi_min,xi_max
   else if (keyword(1:20) .eq. 'PMF_BINS ') then
      read(record,*,iostat=readstat) names,nbins
   else if (keyword(1:20) .eq. 'PMF_METHOD ') then
      read(record,*,iostat=readstat) names,pmf_method
      call upcase(pmf_method)
   else if (keyword(1:20) .eq. 'TS_STRUC ') then
      read(record,*,iostat=readstat) names,ts_file
   else if (keyword(1:20) .eq. 'RECROSS_EQUI ') then
      read(record,*,iostat=readstat) names,recr_equi
   else if (keyword(1:20) .eq. 'CHILD_TOTAL ') then
      read(record,*,iostat=readstat) names,child_tot
   else if (keyword(1:20) .eq. 'CHILD_INTERVAL ') then
      read(record,*,iostat=readstat) names,child_interv
   else if (keyword(1:20) .eq. 'CHILD_PERPOINT ') then
      read(record,*,iostat=readstat) names,child_point
   else if (keyword(1:20) .eq. 'CHILD_EVOL ') then
      read(record,*,iostat=readstat) names,child_evol
   else if (keyword(1:20) .eq. 'RECROSS_MPI ') then
      recross_mpi = .true.
   else if (keyword(1:16) .eq. 'PMF_MINLOC ') then
      read(record,*,iostat=readstat) names,pmf_minloc
   else if (keyword(1:16) .eq. 'PRINT_POLYMER ') then
      read(record,*,iostat=readstat) names,print_poly
   else if (keyword(1:16) .eq. 'RPMD_EN_TOL ') then
      read(record,*,iostat=readstat) names,energy_tol
   else if (keyword(1:16) .eq. 'RECROSS_NOCHECK ') then
      recross_check = .false.
   else if (keyword(1:16) .eq. 'GEN_TEST ') then
      read(record,*,iostat=readstat) names,gen_pr_frac
      gen_pr_act=0
      gen_test=.true.
   else if (keyword(1:16) .eq. 'DEBUG ') then
      do_debug=.true.
      debug_unit=112
   end if
   if (readstat .ne. 0) then
      write(*,*) "ERROR! The qmdff.key file has an invalid formate!"
      write(*,*) "Did you forget keyword-parameters?"
      call fatal
   end if
end do
!
!      Check if all read in parameters are valid.
!      Abort the program if something is not valid
!
if (rank .eq. 0) then
!
!      General sampling trajectory settings
!
   if (gen_step .eq. 0) then
      write(*,*) "No number of start structure generation MD steps defined!"
      write(*,*) "Add the keyword GENERATE_STEPS!"
      call fatal
   else if (equi_step .eq. 0) then
      write(*,*) "No number of umbrella equilibration MD steps defined!"
      write(*,*) "Add the keyword EQUILIBR_STEPS!"
      call fatal
   else if (umbr_step .eq. 0) then
      write(*,*) "No number of umbrella sampling MD steps defined!"
      write(*,*) "Add the keyword UMBRELLA_STEPS!"
      call fatal
   else if (umbr_traj .eq. 0) then
      write(*,*) "No number of single umbrella sampling trajectories per window defined!"
      write(*,*) "Add the keyword UMBRELLA_TRAJS!"
      call fatal
   end if
!
!      Thermostat settings
!
   if (dt .eq. -1.0d0) then
      write(*,*) "No length of time step defined! Add the keyword DELTAT!"
      call fatal
   else if (kelvin .le. 0.0d0) then
      write(*,*) "No temperature defined! Add the keyword TEMP!"
      call fatal
   else if ((thermo .ne. "ANDERSEN") .and. (thermo .ne. "GLE")) then
      write(*,*) "No availiable thermostat was chosen! The Andersen"
      write(*,*) "thermostat as well as the Generalized Langevin Equation (GLE)"
      write(*,*) "thermostat are avaibliable at the moment. Add the keyword THERMOSTAT!"
      call fatal
   end if
!
!      The RPMD beads
!
   if (nbeads .lt. 1) then
      write(*,*) "No valid bead number given! Add the keyword BEAD_NUMBER!"
      call fatal
   end if
!
!     the number of equivalent reaction paths (multiplicator for k(T) value)
!
   if (npaths .lt. 1) then
      write(*,*) "No valid number of equivalent reaction paths given! At least NPATHS 1!"
      call fatal
   end if
!
!     specific umbrella sampling settings
!
   if (k_force .eq. -1.0) then
      write(*,*) "No umbrella force constant defined! Add the keyword UMBRELLA_BIAS!"
      call fatal
   else if ((umbr_lo .eq. -5.0) .or. (umbr_hi .eq. -5.0)) then
      write(*,*) "No range for umbrella samplings defined! Add the keyword UMBRELLA_BONDS!"
      call fatal
   else if (umbr_dist .eq. -1.0) then
      write(*,*) "No distance between umbrella samplings defined! Add the keyword UMBRELLA_DIST!"
      call fatal
   else if (umbr_type .eq. "none") then
      write(*,*) "No type of biased coordinate system is given! Add the keyword UMBRELLA_TYPE!"
      write(*,*) "Coordinate systems currently supported:"
!      write(*,*) "UNIMOLECULAR: A unimolecular reaction: one bond is built, one broken (one molecule)"
      write(*,*) "BIMOLEC: A bimolecular reaction: one bond is built, the other broken."
      write(*,*) "CYCLOADD: A cycloaddition: two bonds are built simultaneously."
      write(*,*) "ADDITION: A bimolecular addition: two bonds are built, one is broken."
      write(*,*) "ADDITION3: A trimolecular addition: three bonds are built, two are broken."
      write(*,*) "ADD3_SOLV: A trimolecular addition: three bonds are built, two are broken, two"
      write(*,*) "           educts are clustered before the reaction (bimol. k(T))."
      write(*,*) "ADDITION4: A tetramolecular addition: four bonds are built, three are broken."
      write(*,*) "ADD4_SOLV: A tetramolecular addition: four bonds are built, three are broken, three"
      write(*,*) "           educts are clustered before the reaction (bimol. k(T))."
      write(*,*) "ATOM_SHIFT: A single atom is transported a certain distance, e.g. on a metal surface."
      write(*,*) "CYCLOREV: A cycloreversion: two bonds are broken simultaneously."
      call fatal
   end if
!
!      prevent usage of unimolecular reactions
!
   if (umbr_type .eq. "UNIMOLECULAR") then
      write(*,*) "Currently, this type is not implemented!"
      call fatal
   end if
!
!      The umbrella integration settings
!
   if (xi_min .eq. 1.0 .or. xi_max .eq. -1.0) then
      write(*,*) "No range for PMF calculation defined! Add the keyword PMF_XI_RANGE!"
      call fatal
   else if (nbins .eq. 0) then
      write(*,*) "No number of windows for PMF calculation defined!"
      write(*,*) "Add the keyword PMF_BINS"
      call fatal
   end if
   call upcase(pmf_minloc)
   write(*,*) "Minloc", pmf_minloc
   if (pmf_minloc .ne. 'ZERO' .and. pmf_minloc .ne. 'PMF_MIN') then
      write(*,*) "The keyword PMF_MINLOC was found but no valid option was"
      write(*,*) "chosen! Take ZERO or PMF_MIN!"
      call fatal
   end if
!
!      The asymptotic educts distance (not needed for ATOM_SHIFT)
!
   if (R_inf .eq. 0.d0) then
     if ((umbr_type .ne. "ATOM_SHIFT") .and. (umbr_type .ne. "CYCLOREV")) then
        write(*,*) "No educt max distance is given! Add the keyword R_INF!"
        call fatal
     end if
   end if
!
!      If no valid TS structure filename was given
!
   if (ts_file .eq. " ") then
      write(*,*) "No filename with the TS structure was given!"
      write(*,*) "Add the keyword TS_STRUC!"
      call fatal
   end if
!
!     check parameters for recrossing calculation
!
   if (recr_equi .eq. 0) then
      write(*,*) "No timestep number for the equilibration of the recrossing parent"
      write(*,*) "trajectory is given! Add the keyword RECROSS_EQUI!"
      call fatal
   else if (child_tot .eq. 0) then
      write(*,*) "No total number of child trajectories for the recrossing calculation"
      write(*,*) "is given! Add the keywod CHILD_TOTAL!"
      call fatal
   else if (child_interv .eq. 0) then
      write(*,*) "No interval between a new couple of child trajectories shall be spawned"
      write(*,*) "is given! Add the keyword CHILD_INTERVAL!"
      call fatal
   else if (child_point .eq. 0) then
      write(*,*) "No number of child trajectories that shall be spawned together at each"
      write(*,*) "spawining time is given! Add the keyword CHILD_PERPOINT!"
      call fatal
   else if (child_evol .eq. 0) then
      write(*,*) "No timestep number for each single child trajectory is given!"
      write(*,*) "Add the keyword CHILD_EVOL!"
   end if
   child_times=child_tot/child_point
   if (child_times .ne. real(child_tot/child_point)) then
      write(*,*) "The total number of child trajectories is no multiplicity of the child"
      write(*,*) "number per sampling point! Alter keywords CHILD_TOTAL and/or CHILD_PERPOINT!"
      call fatal
   end if
!
!      Check print polymer setting for printout of one trajectory: needs to be
!       inside the valid sampling area
!
   if (print_poly .gt. -5.d0) then
      if (print_poly .lt. xi_min .or. print_poly .gt. xi_max) then
         write(*,*) "The 'PRINT_POLYMER' option was activated such that a single umbrella"
         write(*,*) "shall be written out but the Xi location is not in the sampled area!"
         call fatal
      end if
   end if
end if
!
!      Read in dependent information (detailed subinformation of the one given
!      above) or calculate dependent parameters
!
!
!     Define type of thermostat, if the GLE thermostat is used, read in the
!     A (and C) matrix
!
if (thermo .eq. "ANDERSEN") then
   thermostat=0
!
!     Determine the rate at which the thermostat random hits shall be applied:
!     after each andersen_step MD steps (default: 80)
!
   andersen_step=80
   do i = 1, nkey
      next = 1
      record = keyline(i)
      call gettext (record,keyword,next)
      call upcase (keyword)
      string = record(next:120)
      if (keyword(1:20) .eq. 'ANDERSEN_STEP') then
         read(record,*) names,andersen_step
      end if
   end do
else if (thermo .eq. "GLE") then
   if (rank .eq. 0) then
      write(*,*) "The Generalized Langevin Equation (GLE) thermostat is used."
      write(*,*) "Read in the gle_A.txt file with the parameter matrix A..."
   end if
   inquire(file="gle_A.txt",exist=exist)
   if (.not. exist) then
      if (rank .eq. 0) then
         write(*,*) "ERROR! The GLE thermostat was chosen but no gle_A.txt file provided!"
      end if
      call fatal
   end if
   allocate(gle_mat_A(5,5))  ! currently, only 5x5 GLE matrices are needed
   allocate(gle_mat_C(5,5))  ! the GLE C matrix will be set to zero first
   gle_ns=5-1   ! number of dimensions n for GLE thermostat
   open(unit=45,file="gle_A.txt",status="old")
   do i=1,5
      read(45,*,iostat=readstat) gle_mat_A(i,:)
      if (readstat .ne. 0 .and. rank .eq. 0) then
         write(*,*) "ERROR! The GLE matrix file gle_A.txt seems to have a wrong formate!"
         call fatal
      end if
   end do
!
!    convert the GLE matrix entries to atomic units of inverse time
!
   gle_mat_A=gle_mat_A*2.418884326505e-17
   gle_mat_C=0.d0
   close(45)
   thermostat=1
end if

!
!     Option for fixing of distinct atoms in order to make e.g. metal slab calculations
!     possible
!
fix_atoms=.false.
do i = 1, nkey
    next = 1
    record = keyline(i)
    call gettext (record,keyword,next)
    call upcase (keyword)
    string = record(next:120)
    if (keyword(1:20) .eq. 'FIX_ATOMS ') then
       fix_atoms=.true.
       read(record,*) names,fix_file
    end if
    if (fix_atoms) then
       inquire(file=fix_file,exist=exist)
       if (.not. exist) then
          if (rank .eq. 0) then
             write(*,*) "ERROR! The file ",fix_file," with the fixed atoms is not present!"
          end if
          call fatal
       end if
       allocate(fix_list(1000))
!
!     Read in list with fixed atoms
!
       open(unit=213,file=fix_file,status="old")
       k=1
       do
          read(213,*,iostat=readstat) fix_list(k)
          if (readstat .ne. 0) exit
          k=k+1
       end do
       fix_num=k-1
       exit
    end if
end do

!
!    Read in the type of reaction and define the number of forming and breaking bonds
!
call upcase(umbr_type)
if (umbr_type .eq. "UNIMOLECULAR") then
   form_num=1
   break_num=1
   sum_eds=1
else if (umbr_type .eq. "BIMOLEC") then
   form_num=1
   break_num=1
   sum_eds=2
else if (umbr_type .eq. "CYCLOADD") then
   form_num=2
   break_num=0
   sum_eds=2
else if (umbr_type .eq. "MERGING") then
   form_num=1
   break_num=0
   sum_eds=2
else if (umbr_type .eq. "ADDITION") then
   form_num=2
   break_num=1
   sum_eds=2
else if (umbr_type .eq. "ADDITION3") then
   form_num=3
   break_num=2
   sum_eds=3
else if (umbr_type .eq. "ADD3_SOLV") then
   form_num=3
   break_num=2
   sum_eds=2
else if (umbr_type .eq. "ADDITION4") then
   form_num=4
   break_num=3
   sum_eds=4
else if (umbr_type .eq. "ADD4_SOLV") then
   form_num=4
   break_num=3
   sum_eds=2
else if (umbr_type .eq. "ELIMINATION") then
   form_num=1
   break_num=2
   sum_eds=1
else if (umbr_type .eq. "CYCLOREVER") then
   form_num=0
   break_num=2
   sum_eds=1
else if (umbr_type .eq. "REARRANGE") then
   form_num=1
   break_num=1
   sum_eds=1
else if (umbr_type .eq. "DECOM_1BOND") then
   form_num=0
   break_num=1
   sum_eds=1
!
!     If the special ATOM_SHIFT type is used, read in all needed info and go
!     directly to the end of this subsection
!
else if (umbr_type .eq. "ATOM_SHIFT") then
   do i = 1, nkey
      next = 1
      record = keyline(i)
      call gettext (record,keyword,next)
      call upcase (keyword)
      string = record(next:120)
      if (keyword(1:20) .eq. 'SHIFT_ATOM') then
         read(record,*) names,shift_atom
         if (shift_atom .lt. 1 .or. shift_atom .gt. natoms) then
            write(*,*) "The atom listed in SHIFT_ATOM is not in the system!"
         end if
      else if (keyword(1:20) .eq. 'SHIFT_COORD') then
         read(record,*) names,names2
         if (trim(names2) .eq. "X") then
            shift_coord=1
         else if (trim(names2) .eq. "Y") then
            shift_coord=2
         else if (trim(names2) .eq. "Z") then
            shift_coord=3
         else
            write(*,*) "No valid coordinate (x, y or z) is given for SHIFT_COORD!"
            call fatal
         end if
      else if (keyword(1:20) .eq. 'SHIFT_INTERV') then
         read(record,*) names,shift_lo,shift_hi
!
!     convert coordinate borders to bohr
!
         shift_lo=shift_lo/bohr
         shift_hi=shift_hi/bohr
      end if
   end do

   goto 22
else
   if (rank .eq. 0) then
      write(*,*) "This umbrella coordinate set is currently not supported!"
      write(*,*) "Use one of the supported sets (look into the manual)"
   end if
   call fatal
end if

!
!     For unimoleculae reactions: read in structure of educt molecule!
!
do i = 1, nkey
    next = 1
    record = keyline(i)
    call gettext (record,keyword,next)
    call upcase (keyword)
    string = record(next:120)
    if (keyword(1:20) .eq. 'EDUCTS_STRUC ') then
       fix_atoms=.true.
       read(record,*) names,educts_file

       open(unit=31,file=educts_file,status="old",iostat=readstat)
       if (readstat .ne. 0) then
          if (rank .eq. 0) then
             write(*,*) "The file ",educts_file," with the educts structure could not been found!"
             call fatal
          end if
       end if
       allocate(ed_ref(3,natoms))
       read(31,*) ; read(31,*)
       do j=1,natoms
          read(31,*) names,ed_ref(:,j)
       end do
       close(31)
!
!     convert educts structure to bohr
!
       ed_ref=ed_ref/bohr
    end if

end do
!
!     Now read in more detailed information about the possible reaction
!     coordinates and their settings
!
!
!     if the bimolec type is chosen, read in the R_inf parameter as well as atoms of the
!     educts are read in
!     if three or four educts are needed, read them in, too
!
allocate(at_ed(sum_eds,200),n_ed(sum_eds),mass_ed(sum_eds))
at_ed=0
n_ed=0
do i = 1, nkey
   next = 1
   record = keyline(i)
   call gettext (record,keyword,next)
   call upcase (keyword)
   string = record(next:140)
   if (keyword(1:20) .eq. 'EDUCT1 ') then
      read(record,*) names
         k=1
!
!     Read in the atom indices of each single educt: quite complicated in fortran..
!
      act_number=" "
      do j=1,120
         at_index(j)=record(j:j)
         if (at_index(j) .ne. " ") then
            if (len(trim(act_number)) .lt. 1) then
               act_number=at_index(j)
            else
               act_number=trim(act_number) // at_index(j)
            end if
         else
            if (len(trim(act_number)) .ge. 1) then
               read(act_number,*,iostat=istat) at_ed(1,k)
               if (istat .eq. 0) then
                  k=k+1
               else
                  at_ed(1,k)=0
               end if
            end if
            act_number=" "
         end if
      end do
   else if (keyword(1:20) .eq. 'EDUCT2 ') then
      read(record,*) names
      k=1
!
!     Do the same for the second educt..
!
      act_number=" "
      do j=1,120
         at_index(j)=record(j:j)
         if (at_index(j) .ne. " ") then
            if (len(trim(act_number)) .lt. 1) then
               act_number=at_index(j)
            else
               act_number=trim(act_number) // at_index(j)
            end if
         else
            if (len(trim(act_number)) .ge. 1) then
               read(act_number,*,iostat=istat) at_ed(2,k)
               if (istat .eq. 0) then
                  k=k+1
               else
                  at_ed(2,k)=0
               end if
            end if
            act_number=" "
         end if
      end do
   else if (keyword(1:20) .eq. 'EDUCT3 ') then
      if (sum_eds .gt. 2) then
         read(record,*) names
         k=1
!
!     Do the same for the third educt..
!
         act_number=" "
         do j=1,120
            at_index(j)=record(j:j)
            if (at_index(j) .ne. " ") then
               if (len(trim(act_number)) .lt. 1) then
                  act_number=at_index(j)
               else
                  act_number=trim(act_number) // at_index(j)
               end if
            else
               if (len(trim(act_number)) .ge. 1) then
                  read(act_number,*,iostat=istat) at_ed(3,k)
                  if (istat .eq. 0) then
                     k=k+1
                  else
                     at_ed(3,k)=0
                  end if
               end if
               act_number=" "
            end if
         end do
      end if
   else if (keyword(1:20) .eq. 'EDUCT4 ') then
      if (sum_eds .eq. 4) then
         read(record,*) names
         k=1
!
!     Do the same for the fourth educt..
!
         act_number=" "
         do j=1,120
            at_index(j)=record(j:j)
            if (at_index(j) .ne. " ") then
               if (len(trim(act_number)) .lt. 1) then
                  act_number=at_index(j)
               else
                  act_number=trim(act_number) // at_index(j)
               end if
            else
               if (len(trim(act_number)) .ge. 1) then
                  read(act_number,*,iostat=istat) at_ed(4,k)
                  if (istat .eq. 0) then
                     k=k+1
                  else
                     at_ed(4,k)=0
                  end if
               end if
               act_number=" "
            end if
         end do
      end if
   end if
end do
!
!     Control if additional settings are correct/useful/appropriate
!
do i=1,sum_eds
   if (at_ed(i,1) .eq. 0) then
      if (rank .eq. 0) then
         write(*,'(a,i1,a,i1,a)') "No atom for educt ",i," is given! &
              &  Add the keyword EDUCT",i,"!"
         call fatal
      end if
   end if
end do
!
!     Control the elements of the educts arrays for doublings, emptyness etc..
!     for all educts in one loop!
!
do l=1,sum_eds
   active=.true.
   i=0
   do while (active)
      i=i+1
      if (at_ed(l,i) .eq. 0) then
         active=.false.
      else if ((at_ed(l,i) .gt. natoms) .or. (at_ed(l,i) .lt. 0)) then
         if (rank .eq. 0) then
            write(*,'(a,i1,a)') "Educt ",l," has an atom out of range!"
            call fatal
         end if
      else
         do p=1,l
            do j=1,n_ed(p)
               if (at_ed(l,i) .eq. at_ed(p,j)) then
                  if (rank .eq. 0) then
                    write(*,'(a,i4,a,i1,a,i1,a)') " The atom No. ",at_ed(l,i), &
                        & " is shared between educt ",l," and ",p,"!"
                    call fatal
                  end if
               end if
            end do
         end do
      end if
   end do
   if (i .eq. 1) then
      if (rank .eq. 0) then
         write(*,*) "One of the educts has no atoms!"
         call fatal
      end if
   end if
   n_ed(l)=i-1
end do
if (sum(n_ed) .ne. natoms) then
   if (rank .eq. 0) then
      write(*,*) "The sum of educt atoms isn´t equal to the total number of atoms!"
      call fatal
   end if
end if
!     allocate arrays for breaking and forming bonds as well as for its
!     reading in
!

allocate(bond_form(form_num,2))
allocate(bond_break(break_num,2))
allocate(input_form(form_num))
allocate(input_break(break_num))
!
!     Read in the forming and breaking bonds to define the reaction coordinate Xi
!
do i = 1, nkey
   next = 1
   record = keyline(i)
   call gettext (record,keyword,next)
   call upcase (keyword)
   string = record(next:120)
   if (keyword(1:20) .eq. 'BOND_FORM ') then
!
!     Read all bonds with one command
!
      read(record,*) names,input_form
      do j=1,form_num
         l=1
         act_bond=input_form(j)
!
!     Analyze each substring, extract the numbers of the two participating atoms
!
         act_number=" "
         do k=1,20
            act_char=act_bond(k:k)
            if (act_char .ne. "-" .and. act_char .ne. " ") then
               act_number=trim(act_number) // act_char
            else
               if (len(trim(act_number)) .ne. 0) then
                  read(act_number,*) bond_form(j,l)
                  l=l+1
                  act_number=" "
               end if
            end if
!
!     if both bonded parters are found, exit the inner loop
!
            if (l .gt. 2) exit
         end do
      end do
   else if (keyword(1:20) .eq. 'BOND_BREAK ') then
      read(record,*) names,input_break
      do j=1,break_num
         l=1
         act_bond=input_break(j)
!
!     Analyze each substring, extract the numbers of the two participating atoms
!
         act_number=" "
         do k=1,20
            act_char=act_bond(k:k)
            if (act_char .ne. "-" .and. act_char .ne. " ") then
               act_number=trim(act_number) // act_char
            else
               if (len(trim(act_number)) .ne. 0) then
                  read(act_number,*) bond_break(j,l)
                  l=l+1
                  act_number=" "
               end if
            end if
!
!     if both bonded parters are found, exit the inner loop
!
            if (l .gt. 2) exit
         end do
      end do
   end if
end do

!
!     Read in the forming and breaking bonds to define the reaction coordinate Xi
!
do i = 1, nkey
   next = 1
   record = keyline(i)
   call gettext (record,keyword,next)
   call upcase (keyword)
   string = record(next:120)
   if (keyword(1:20) .eq. 'BOND_FORM ') then
!
!     Read all bonds with one command
!
      read(record,*) names,input_form
      do j=1,form_num
         l=1
         act_bond=input_form(j)
!
!     Analyze each substring, extract the numbers of the two participating atoms
!
         act_number=" "
         do k=1,20
            act_char=act_bond(k:k)
            if (act_char .ne. "-" .and. act_char .ne. " ") then
               act_number=trim(act_number) // act_char
            else
               if (len(trim(act_number)) .ne. 0) then
                  read(act_number,*) bond_form(j,l)
                  l=l+1
                  act_number=" "
               end if
            end if
!
!     if both bonded parters are found, exit the inner loop
!
            if (l .gt. 2) exit
         end do
      end do
   else if (keyword(1:20) .eq. 'BOND_BREAK ') then
      read(record,*) names,input_break
      do j=1,break_num
         l=1
         act_bond=input_break(j)
!
!     Analyze each substring, extract the numbers of the two participating atoms
!
         act_number=" "
         do k=1,20
            act_char=act_bond(k:k)
            if (act_char .ne. "-" .and. act_char .ne. " ") then
               act_number=trim(act_number) // act_char
            else
               if (len(trim(act_number)) .ne. 0) then
                  read(act_number,*) bond_break(j,l)
                  l=l+1
                  act_number=" "
               end if
            end if
!
!     if both bonded parters are found, exit the inner loop
!
            if (l .gt. 2) exit
         end do
      end do
   end if
end do
!
!      Continue here for the ATOM_SHIFT coordinate
!
22 continue
!
!      Check methods for PMF integration. If none is set, set it to integration
!        as default
!
call upcase(pmf_method)
if ((pmf_method .ne. 'WHAM') .and. (pmf_method .ne. 'INTEGRATION') .and. &
              & (pmf_method .ne. 'ALL')) then
   if (rank .eq. 0) then
      write(15,*) "No (existent) method for PMF calculation is given by the user!"
      write(15,*) "Umbrella integration will be used instead."
   end if
   pmf_method="INTEGRATION"
end if

!
!      Read in the frequency of structure printouts for debug trajectory
!      If no keyword given, take 1 fs interval as default
!
if (print_poly .ne. -5.d0) then
   do i = 1, nkey
      next = 1
      record = keyline(i)
      call gettext (record,keyword,next)
      call upcase (keyword)
      string = record(next:120)
      if (keyword(1:11) .eq. 'TDUMP ') then
         read(record,*) names,dtdump
!
!     calculate the frequency of print outs
!
         iwrite = nint(dtdump/dt_info*0.001d0)
      end if
   end do
end if

!
!     For MECHANOCHEMISTRY simulations: Add forces to certain atoms!
!
!     Main keyword: ADD_FORCE
!
add_force = .false.
do i = 1, nkey
    next = 1
    record = keyline(i)
    call gettext (record,keyword,next)
    call upcase (keyword)
    string = record(next:120)
    if (keyword(1:11) .eq. 'ADD_FORCE ') then
       add_force=.true.
       write(*,*) "The keyword ADD_FORCE was found!"
       write(*,*) "Therefore a mechanochemistry trajectory will be"
       write(*,*) "calculated, until a rupture event is monitored and"
       write(*,*) "the user will be informed."
       write(*,*) "The keywords FORCE1 and FORCE2 must be present, which the formate:"
       write(*,*) "FORCE1/2 [atom index] [force constant (N)] [vector (x,y,z)]"
       exit
    end if
end do

!
!     If Mechanochemistry was activated: read in the atoms as well as
!      vectors and strenghts for the applied forces
!
if (add_force) then
   do i = 1, nkey
      next = 1
      record = keyline(i)
      call gettext (record,keyword,next)
      call upcase (keyword)
      string = record(next:120)
!
!     Ordering: [Atom index] [Force (N)] [Vector (x,y,z)]
!
      if (keyword(1:11) .eq. 'FORCE1 ') then
         read(record,*) names,force1_at,force1_k,force1_v(1),force1_v(2),force1_v(3)
      else if (keyword(1:11) .eq. 'FORCE2 ') then
         read(record,*) names,force2_at,force2_k,force2_v(1),force2_v(2),force2_v(3)
      end if
   end do
!
!     normalize the force vectors
!
   force1_v=force1_v/(sqrt(force1_v(1)**2+force1_v(2)**2+force1_v(3)*2))
   force2_v=force2_v/(sqrt(force2_v(1)**2+force2_v(2)**2+force2_v(3)*2))
!
!     transform units of force constants from Newton to atomic units
!
   if ((force1_at .eq. 0) .or. (force1_at .gt. natoms) .or.&
              &  (force2_at .eq. 0) .or. (force2_at .gt. natoms)) then
      write(*,*) "ERROR! The atoms at which the force shall be applied are"
      write(*,*) " not part of the system!"
      call fatal
   end if
  !  write(*,*) "FORCE setup:"
  !  write(*,*) "At atom",force1_at," a force of",force1_k," N is applied."
  !  write(*,*) "At atom",force2_at," a force of",force2_k," N is applied."
   write(*,*) newton2au
   force1_k=force1_k*newton2au
   force2_k=force2_k*newton2au
!
!     Open file in which distances between the atoms with attached force will
!     be written
!
   write(*,*) "Distances between atoms with attached force will be written to"
   write(*,*) "file force_distances.dat"
   open(unit=78,file="force_distances.dat",status="replace")
   write(78,*) "# distances between the atoms with attached force, in Angstrom:"
end if

end subroutine rpmd_read

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: chrinide/EVB-QMDFF

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: chrinide/EVB-QMDFF

max_forks_repo_licenses: MIT Unlicense


text:
Candice did hair in Davis from 2007 to 2012. She is now working at Jane Salon in Sacramento. Her passions are fun cuts and creative color for men and women. She also specializes in fusion extensions, which are individual pieces of 100% human hair fused to your natural hair with a keratin protein bond that lasts 36 months.

Her Prices Were
Haircuts starting @ $35
Color starting @ $55
Highlight starting @ $85
Blowout $35
Deep condition $20
Deep cond/blowout combo $40
Fusion extensions price upon request

20110713 15:44:08 nbsp I have seen Candace twice now for a haircut and she does excellent work!  She is very experienced with classic and trendy haircuts.  Is it easy to get an appointment with her and she charges some of the most reasonable rates in Davis. Users/JoyRobbins

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search