hdiff output

r30609/commons.f90 2016-07-06 15:36:10.335694852 +0100 r30608/commons.f90 2016-07-06 15:36:12.171719678 +0100
110:      &        HARMONICDONTMOVE, DUMPUNIQUE, FREEZESAVE, TBP, RBSYMT, PTMCDUMPSTRUCT, PTMCDUMPENERT, PYCOLDFUSION, MONITORT,&110:      &        HARMONICDONTMOVE, DUMPUNIQUE, FREEZESAVE, TBP, RBSYMT, PTMCDUMPSTRUCT, PTMCDUMPENERT, PYCOLDFUSION, MONITORT,&
111:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, SANDBOXT, &111:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, SANDBOXT, &
112:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &112:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &
113:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &113:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &
114:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &114:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &
115:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &115:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &
116:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, & 116:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, & 
117:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &117:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &
118:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &118:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &
119:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &119:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &
120:      &        MLP3T, MKTRAPT, MLPB3T, MLPB3NEWT, MULTIPOTT, QCIAMBERT, MLPNEWREG, DJWRBT, STEALTHYT, LJADDT, QCINOREPINT, RIGIDMDT120:      &        MLP3T, MKTRAPT, MLPB3T, MLPB3NEWT, MULTIPOTT, QCIAMBERT, MLPNEWREG, DJWRBT, STEALTHYT, LJADDT, QCINOREPINT
121: !121: !
122:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:) 122:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:) 
123:       DOUBLE PRECISION, ALLOCATABLE:: ATMASS(:)123:       DOUBLE PRECISION, ALLOCATABLE:: ATMASS(:)
124:       DOUBLE PRECISION, ALLOCATABLE:: SPECMASS(:) 124:       DOUBLE PRECISION, ALLOCATABLE:: SPECMASS(:) 
125: 125: 
126: ! csw34> FREEZEGROUP variables126: ! csw34> FREEZEGROUP variables
127: !127: !
128:       INTEGER :: GROUPCENTRE128:       INTEGER :: GROUPCENTRE
129:       DOUBLE PRECISION :: GROUPRADIUS129:       DOUBLE PRECISION :: GROUPRADIUS
130:       CHARACTER (LEN=2) :: FREEZEGROUPTYPE130:       CHARACTER (LEN=2) :: FREEZEGROUPTYPE


r30609/keywords.f 2016-07-06 15:36:10.695699721 +0100 r30608/keywords.f 2016-07-06 15:36:12.551724814 +0100
794:       COULN=0794:       COULN=0
795:       COULQ=0.0D0795:       COULQ=0.0D0
796:       COULSWAP = 0.0D0796:       COULSWAP = 0.0D0
797:       COULTEMP = 0.0D0797:       COULTEMP = 0.0D0
798:       GAYBERNET=.FALSE.798:       GAYBERNET=.FALSE.
799:       ELLIPSOIDT=.FALSE.799:       ELLIPSOIDT=.FALSE.
800:       PYGPERIODICT=.FALSE.800:       PYGPERIODICT=.FALSE.
801:       LJCAPSIDT=.FALSE.801:       LJCAPSIDT=.FALSE.
802:       PYBINARYT=.FALSE.802:       PYBINARYT=.FALSE.
803:       pyt = .false.803:       pyt = .false.
804:       RIGIDMDT = .false. 
805:       PYOVERLAPTHRESH=1.0D0804:       PYOVERLAPTHRESH=1.0D0
806:       PYCFTHRESH=0.1D0 ! cold fusion threshold805:       PYCFTHRESH=0.1D0 ! cold fusion threshold
807:       LJSITE=.FALSE.806:       LJSITE=.FALSE.
808:       SWAPMOVEST=.FALSE.807:       SWAPMOVEST=.FALSE.
809:       STICKYT=.FALSE.808:       STICKYT=.FALSE.
810:       RIGID=.FALSE.809:       RIGID=.FALSE.
811:       TIPID=4810:       TIPID=4
812:       HEIGHT=0.5D0811:       HEIGHT=0.5D0
813:       OTPT=.FALSE.812:       OTPT=.FALSE.
814:       LJMFT=.FALSE.813:       LJMFT=.FALSE.
6511:          CALL READF(GBNU)6510:          CALL READF(GBNU)
6512:          CALL READF(SIGNOT)6511:          CALL READF(SIGNOT)
6513:          CALL READF(EPSNOT)6512:          CALL READF(EPSNOT)
6514: !         ALLOCATE(SITE(NRBSITES,3))6513: !         ALLOCATE(SITE(NRBSITES,3))
6515: 6514: 
6516:          LPL    = 0.5D0 * SIGNOT * GBKAPPA6515:          LPL    = 0.5D0 * SIGNOT * GBKAPPA
6517:          LPR    = 0.5D0 * SIGNOT6516:          LPR    = 0.5D0 * SIGNOT
6518:          LPRSQ  = LPR * LPR6517:          LPRSQ  = LPR * LPR
6519:          LSQDFR = LPL * LPL - LPRSQ6518:          LSQDFR = LPL * LPL - LPRSQ
6520: 6519: 
6521:       ELSE IF (WORD .EQ. 'RIGIDMD') THEN 
6522:          ! Syntax:  
6523:  
6524:          RIGIDMDT = .TRUE. 
6525:   
6526:       ELSE IF (WORD .EQ. 'PY') THEN6520:       ELSE IF (WORD .EQ. 'PY') THEN
6527:          ! Syntax: PY sig_0 eps_0 [cut] [XYZ boxx boxy boxz]6521:          ! Syntax: PY sig_0 eps_0 [cut] [XYZ boxx boxy boxz]
6528: 6522: 
6529:          PYT = .TRUE.6523:          PYT = .TRUE.
6530:          CALL READF(PYSIGNOT)6524:          CALL READF(PYSIGNOT)
6531:          CALL READF(PYEPSNOT)6525:          CALL READF(PYEPSNOT)
6532: 6526: 
6533:          RIGID = .TRUE.6527:          RIGID = .TRUE.
6534:          ! Rigid body SITE, NRBSITES, NTSITES information specified in py_input routine6528:          ! Rigid body SITE, NRBSITES, NTSITES information specified in py_input routine
6535: 6529: 


r30609/mc.F 2016-07-06 15:36:11.075704859 +0100 r30608/mc.F 2016-07-06 15:36:12.903729575 +0100
962:                 CALL TAKESTEPGB(JP)962:                 CALL TAKESTEPGB(JP)
963: 963: 
964:                ELSE IF (PYT) THEN964:                ELSE IF (PYT) THEN
965: !               seed the random number generator with system time + MYNODE (for MPI runs)965: !               seed the random number generator with system time + MYNODE (for MPI runs)
966:                   IF(RANDOMSEEDT) THEN966:                   IF(RANDOMSEEDT) THEN
967:                      CALL DATE_AND_TIME(datechar,timechar,zonechar,values)967:                      CALL DATE_AND_TIME(datechar,timechar,zonechar,values)
968:                      itime1= values(6)*60 + values(7)968:                      itime1= values(6)*60 + values(7)
969:                      CALL SDPRND(itime1+MYNODE)969:                      CALL SDPRND(itime1+MYNODE)
970:                   END IF970:                   END IF
971:                 call py_takestep(jp)971:                 call py_takestep(jp)
972:                 if (rigidmdt) call rigidmd_nvt 
973: 972: 
974:                ELSE IF (PYGPERIODICT.OR.PYBINARYT.OR.LJCAPSIDT) THEN973:                ELSE IF (PYGPERIODICT.OR.PYBINARYT.OR.LJCAPSIDT) THEN
975: !               seed the random number generator with system time + MYNODE (for MPI runs)974: !               seed the random number generator with system time + MYNODE (for MPI runs)
976:                   IF(RANDOMSEEDT) THEN975:                   IF(RANDOMSEEDT) THEN
977:                      CALL DATE_AND_TIME(datechar,timechar,zonechar,values)976:                      CALL DATE_AND_TIME(datechar,timechar,zonechar,values)
978:                      itime1= values(6)*60 + values(7)977:                      itime1= values(6)*60 + values(7)
979:                      CALL SDPRND(itime1+MYNODE)978:                      CALL SDPRND(itime1+MYNODE)
980:                   END IF979:                   END IF
981:                   IF(SWAPMOVEST) THEN980:                   IF(SWAPMOVEST) THEN
982:                      CALL TAKESTEPSWAPMOVES(JP)981:                      CALL TAKESTEPSWAPMOVES(JP)


r30609/py.f90 2016-07-06 15:36:11.467710165 +0100 r30608/py.f90 2016-07-06 15:36:13.263734444 +0100
 89:             ell%shapemat_bodyframe(k, k) = 1.d0 / (ell%semiaxes(k) ** 2) 89:             ell%shapemat_bodyframe(k, k) = 1.d0 / (ell%semiaxes(k) ** 2)
 90:             ell%shapemat_det = ell%shapemat_det * ell%shapemat_bodyframe(k, k)  90:             ell%shapemat_det = ell%shapemat_det * ell%shapemat_bodyframe(k, k) 
 91:         end do 91:         end do
 92:  92: 
 93:         ell%shapemat_inv_det = 1.d0 / ell%shapemat_det 93:         ell%shapemat_inv_det = 1.d0 / ell%shapemat_det
 94:         call rmdrvt(site%p, site%rotmat, dummy, dummy, dummy, .false.) 94:         call rmdrvt(site%p, site%rotmat, dummy, dummy, dummy, .false.)
 95:         ell%shapemat_molframe(:, :) = matmul(site%rotmat, matmul(ell%shapemat_bodyframe, transpose(site%rotmat))) 95:         ell%shapemat_molframe(:, :) = matmul(site%rotmat, matmul(ell%shapemat_bodyframe, transpose(site%rotmat)))
 96:  96: 
 97:     end subroutine initialize_py_ellipsoid 97:     end subroutine initialize_py_ellipsoid
 98:  98: 
 99:     subroutine pairwise_py(sitei, sitej, grad_contrib, energy_contrib, gtest, t_rep, t_att) 99:     subroutine pairwise_py(sitei, sitej, grad_contrib, energy_contrib, gtest)
100:         ! compute energy and gradient due to a pair of sites100:         ! compute energy and gradient due to a pair of sites
101: 101: 
102:         use commons, only: pysignot, paramonovcutoff, pyepsnot102:         use commons, only: pysignot, paramonovcutoff, pyepsnot
103: 103: 
104:         implicit none104:         implicit none
105:         type(py_site), target, intent(inout) :: sitei, sitej105:         type(py_site), target, intent(inout) :: sitei, sitej
106:         logical, intent(in) :: gtest106:         logical, intent(in) :: gtest
107:         double precision, intent(out) :: energy_contrib, grad_contrib(12)107:         double precision, intent(out) :: energy_contrib, grad_contrib(12)
108: 108: 
109:         type(py_ellipsoid), pointer :: elli, ellj   ! pointers to ellipsoids inside the sites109:         type(py_ellipsoid), pointer :: elli, ellj   ! pointers to ellipsoids inside the sites
110:         double precision :: g, g_pm1, g_pm2, g_pm6  ! powers of g110:         double precision :: g, g_pm1, g_pm2, g_pm6  ! powers of g
111:         double precision :: F_pm1_2, dg_dF, dg_dr(3), dg_drmod  ! F_pm1_2 = F ** (-1/2)111:         double precision :: F_pm1_2, dg_dF, dg_dr(3), dg_drmod  ! F_pm1_2 = F ** (-1/2)
112:         double precision :: e_rep, g_rep(12), e_att, g_att(12)112:         double precision :: e_rep, g_rep(12), e_att, g_att(12)
113:         double precision :: g_pm7, g_pm12, g_pm13, g_pm14113:         double precision :: g_pm7, g_pm12, g_pm13, g_pm14
114:         double precision :: dU_dF, dU_drmod, torquei(6), torquej(6), t_att(6), t_rep(6)114:         double precision :: dU_dF, dU_drmod
115:          
116: 115: 
117:         ! repulsive part116:         ! repulsive part
118:         elli => sitei%rep117:         elli => sitei%rep
119:         ellj => sitej%rep118:         ellj => sitej%rep
120: 119: 
121:         call compute_py_values(elli, ellj, gtest, g, g_pm1, g_pm2, g_pm6, F_pm1_2, dg_dF, dg_dr, dg_drmod, torquei, torquej)120:         call compute_py_values(elli, ellj, gtest, g, g_pm1, g_pm2, g_pm6, F_pm1_2, dg_dF, dg_dr, dg_drmod)
122:         if(paramonovcutoff .and. above_cutoff) then121:         if(paramonovcutoff .and. above_cutoff) then
123:             e_rep = 0.d0122:             e_rep = 0.d0
124:             g_rep(:) = 0.d0123:             g_rep(:) = 0.d0
125:         else ! not over the cutoff124:         else ! not over the cutoff
126:             g_pm12 = g_pm6 * g_pm6125:             g_pm12 = g_pm6 * g_pm6
127:             g_pm13 = g_pm12 * g_pm1126:             g_pm13 = g_pm12 * g_pm1
128: 127: 
129:             e_rep = g_pm12128:             e_rep = g_pm12
130:             if(paramonovcutoff) e_rep = e_rep + gcut_pm12 * (6.d0 * gcut_pm2 / g_pm2 - 7.d0)129:             if(paramonovcutoff) e_rep = e_rep + gcut_pm12 * (6.d0 * gcut_pm2 / g_pm2 - 7.d0)
131: 130: 
138:                 g_rep(10: 12) = dU_dF * ellj%dF_dp(:) + dU_drmod * ellj%dr_dp(:)137:                 g_rep(10: 12) = dU_dF * ellj%dF_dp(:) + dU_drmod * ellj%dr_dp(:)
139: 138: 
140:                 if(paramonovcutoff) then139:                 if(paramonovcutoff) then
141:                     ! add splines to achieve cutoff140:                     ! add splines to achieve cutoff
142:                     g_rep(1: 3) = g_rep(1: 3) + 2.d0 * gcut_pm14 * g * dg_dr(:)141:                     g_rep(1: 3) = g_rep(1: 3) + 2.d0 * gcut_pm14 * g * dg_dr(:)
143:                     g_rep(7: 9) = g_rep(7: 9) + 2.d0 * gcut_pm14 * g * (dg_dF * elli%dF_dp(:) + dg_drmod * elli%dr_dp(:))142:                     g_rep(7: 9) = g_rep(7: 9) + 2.d0 * gcut_pm14 * g * (dg_dF * elli%dF_dp(:) + dg_drmod * elli%dr_dp(:))
144:                     g_rep(10: 12) = g_rep(10: 12) + 2.d0 * gcut_pm14 * g * (dg_dF * ellj%dF_dp(:) + dg_drmod * ellj%dr_dp(:))143:                     g_rep(10: 12) = g_rep(10: 12) + 2.d0 * gcut_pm14 * g * (dg_dF * ellj%dF_dp(:) + dg_drmod * ellj%dr_dp(:))
145:                 end if144:                 end if
146: 145: 
147:                 g_rep(4: 6) = -g_rep(1: 3)  ! (CW eq 34)146:                 g_rep(4: 6) = -g_rep(1: 3)  ! (CW eq 34)
148:                 t_rep(1:3)=24*pyepsnot*sitei%rep%strength*torquei(1:3) 
149:                 t_rep(4:6)=24*pyepsnot*sitej%rep%strength*torquej(1:3) 
150:             end if147:             end if
151:         end if148:         end if
152: 149: 
153:         ! attractive part: if ellipsoids are the same, keep the old values and ellipsoid pointers150:         ! attractive part: if ellipsoids are the same, keep the old values and ellipsoid pointers
154:         if(.not.(sitei%ells_same .and. sitej%ells_same)) then151:         if(.not.(sitei%ells_same .and. sitej%ells_same)) then
155:             elli => sitei%att152:             elli => sitei%att
156:             ellj => sitej%att153:             ellj => sitej%att
157:             call compute_py_values(elli, ellj, gtest, g, g_pm1, g_pm2, g_pm6, F_pm1_2, dg_dF, dg_dr, dg_drmod, torquei, torquej)154:             call compute_py_values(elli, ellj, gtest, g, g_pm1, g_pm2, g_pm6, F_pm1_2, dg_dF, dg_dr, dg_drmod)
158:         end if155:         end if
159: 156: 
160:         if(paramonovcutoff .and. above_cutoff) then157:         if(paramonovcutoff .and. above_cutoff) then
161:             e_att = 0.d0158:             e_att = 0.d0
162:             g_att(:) = 0.d0159:             g_att(:) = 0.d0
163:             160:             
164:         else ! not over cutoff161:         else ! not over cutoff
165:             g_pm12 = g_pm6 * g_pm6162:             g_pm12 = g_pm6 * g_pm6
166:             g_pm13 = g_pm12 * g_pm1163:             g_pm13 = g_pm12 * g_pm1
167: 164: 
177:                 g_att(7: 9) = dU_dF * elli%dF_dp(:) + dU_drmod * elli%dr_dp(:)174:                 g_att(7: 9) = dU_dF * elli%dF_dp(:) + dU_drmod * elli%dr_dp(:)
178:                 g_att(10: 12) = dU_dF * ellj%dF_dp(:) + dU_drmod * ellj%dr_dp(:)175:                 g_att(10: 12) = dU_dF * ellj%dF_dp(:) + dU_drmod * ellj%dr_dp(:)
179: 176: 
180:                 if(paramonovcutoff) then177:                 if(paramonovcutoff) then
181:                     g_att(1: 3) = g_att(1: 3) - gcut_pm8 * g * dg_dr(:)178:                     g_att(1: 3) = g_att(1: 3) - gcut_pm8 * g * dg_dr(:)
182:                     g_att(7: 9) = g_att(7: 9) - gcut_pm8 * g * (dg_dF * elli%dF_dp(:) + dg_drmod * elli%dr_dp(:))179:                     g_att(7: 9) = g_att(7: 9) - gcut_pm8 * g * (dg_dF * elli%dF_dp(:) + dg_drmod * elli%dr_dp(:))
183:                     g_att(10: 12) = g_att(10: 12) - gcut_pm8 * g * (dg_dF * ellj%dF_dp(:) + dg_drmod * ellj%dr_dp(:))180:                     g_att(10: 12) = g_att(10: 12) - gcut_pm8 * g * (dg_dF * ellj%dF_dp(:) + dg_drmod * ellj%dr_dp(:))
184:                 end if181:                 end if
185: 182: 
186:                 g_att(4: 6) = -g_att(1: 3)183:                 g_att(4: 6) = -g_att(1: 3)
187:                 t_att(1:3)=12*pyepsnot*sitei%att%strength*torquei(4:6) 
188:                 t_att(4:6)=12*pyepsnot*sitej%att%strength*torquej(4:6) 
189:             end if184:             end if
190:         end if185:         end if
191: 186: 
192:         ! add on final factors187:         ! add on final factors
193:         energy_contrib = 4.d0 * pyepsnot * (sitei%rep%strength * sitej%rep%strength * e_rep &188:         energy_contrib = 4.d0 * pyepsnot * (sitei%rep%strength * sitej%rep%strength * e_rep &
194:             & + sitei%att%strength * sitej%att%strength * e_att)189:             & + sitei%att%strength * sitej%att%strength * e_att)
195:         if(gtest) grad_contrib(:) = 24.d0 * pyepsnot * (sitei%rep%strength * sitej%rep%strength * g_rep(:) &190:         if(gtest) grad_contrib(:) = 24.d0 * pyepsnot * (sitei%rep%strength * sitej%rep%strength * g_rep(:) &
196:             & + sitei%att%strength * sitej%att%strength * g_att(:))191:             & + sitei%att%strength * sitej%att%strength * g_att(:))
197: 192: 
198:     end subroutine pairwise_py193:     end subroutine pairwise_py
199: 194: 
200:     subroutine compute_py_values(elli, ellj, gtest, g, g_pm1, g_pm2, g_pm6, F_pm1_2, dg_dF, dg_dr, dg_drmod, torquei, torquej)195:     subroutine compute_py_values(elli, ellj, gtest, g, g_pm1, g_pm2, g_pm6, F_pm1_2, dg_dF, dg_dr, dg_drmod)
201:         use commons, only: pysignot, paramonovcutoff196:         use commons, only: pysignot, paramonovcutoff
202:         implicit none197:         implicit none
203: 198: 
204:         type(py_ellipsoid), target, intent(inout) :: elli, ellj199:         type(py_ellipsoid), target, intent(inout) :: elli, ellj
205:         logical, intent(in) :: gtest200:         logical, intent(in) :: gtest
206:         double precision, intent(out) :: g, g_pm1, g_pm2, g_pm6, F_pm1_2, dg_dF, dg_dr(3), dg_drmod201:         double precision, intent(out) :: g, g_pm1, g_pm2, g_pm6, F_pm1_2, dg_dF, dg_dr(3), dg_drmod
207: 202: 
208:         type(py_site), pointer :: sitei, sitej203:         type(py_site), pointer :: sitei, sitej
209:         type(py_molecule), pointer :: moli, molj204:         type(py_molecule), pointer :: moli, molj
210:         double precision :: lambda_c, F, rij(3), rij_p2, rijmod, rij_hat(3)205:         double precision :: lambda_c, F, rij(3), rij_p2, rijmod, rij_hat(3)
211:         double precision :: x_c(3), xc_minus_ri(3), xc_minus_rj(3)  ! PY eq 9; last term in CW eq 28206:         double precision :: x_c(3), xc_minus_ri(3), xc_minus_rj(3)  ! PY eq 9; last term in CW eq 28
212:         double precision :: xlambdac_term1(3, 3), xlambdac_term2(3) ! PY eq 5 with lambda_c from after PY eq 7207:         double precision :: xlambdac_term1(3, 3), xlambdac_term2(3) ! PY eq 5 with lambda_c from after PY eq 7
213:         double precision, intent(out) :: torquei(6), torquej(6) 
214:         integer :: i208:         integer :: i
215: 209: 
216:         sitei => elli%site210:         sitei => elli%site
217:         sitej => ellj%site211:         sitej => ellj%site
218:         moli => sitei%molecule212:         moli => sitei%molecule
219:         molj => sitej%molecule213:         molj => sitej%molecule
220: 214: 
221:         call compute_contact(elli, ellj, lambda_c, F)215:         call compute_contact(elli, ellj, lambda_c, F)
222: 216: 
223:         rij(:) = elli%r(:) - ellj%r(:)217:         rij(:) = elli%r(:) - ellj%r(:)
255: 249: 
256:                 dg_dF = 0.5d0 * rijmod * F_pm1_2 / (F * pysignot)250:                 dg_dF = 0.5d0 * rijmod * F_pm1_2 / (F * pysignot)
257:                 dg_dr(:) = (1.d0 - F_pm1_2) * rij_hat(:) / pysignot + dg_dF * elli%dF_dr(:) ! CW eq 17 first term251:                 dg_dr(:) = (1.d0 - F_pm1_2) * rij_hat(:) / pysignot + dg_dF * elli%dF_dr(:) ! CW eq 17 first term
258: 252: 
259:                 if(paramonovcutoff) dg_drmod = (1.d0 - F_pm1_2) / pysignot253:                 if(paramonovcutoff) dg_drmod = (1.d0 - F_pm1_2) / pysignot
260: 254: 
261:                 do i = 1, 3255:                 do i = 1, 3
262:                     elli%dr_dp(i) = dot_product(rij_hat(:), matmul(moli%rotmat_deriv(i, :, :), sitei%r_molframe(:)))256:                     elli%dr_dp(i) = dot_product(rij_hat(:), matmul(moli%rotmat_deriv(i, :, :), sitei%r_molframe(:)))
263:                     ellj%dr_dp(i) = dot_product(-rij_hat(:), matmul(molj%rotmat_deriv(i, :, :), sitej%r_molframe(:)))257:                     ellj%dr_dp(i) = dot_product(-rij_hat(:), matmul(molj%rotmat_deriv(i, :, :), sitej%r_molframe(:)))
264:                 end do258:                 end do
265: ! torques (PY eq C11) 
266:             torquei(1:3)=-2.d0*g_pm6*g_pm6*g_pm1*dg_dF*(cross(xc_minus_ri,elli%dF_dr)) 
267:             torquei(4:6)=g_pm6*g_pm1*dg_dF*(cross(xc_minus_ri,elli%dF_dr)) 
268:             torquej(1:3)=-2.d0*g_pm6*g_pm6*g_pm1*dg_dF*(cross(xc_minus_rj,ellj%dF_dr)) 
269:             torquej(4:6)=g_pm6*g_pm1*dg_dF*(cross(xc_minus_rj,ellj%dF_dr)) 
270:  
271:             end if259:             end if
272:         end if260:         end if
273: contains 
274: FUNCTION cross(a, b) 
275:   DOUBLE PRECISION, DIMENSION(3) :: cross 
276:   DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: a, b 
277:  
278:   cross(1) = a(2) * b(3) - a(3) * b(2) 
279:   cross(2) = a(3) * b(1) - a(1) * b(3) 
280:   cross(3) = a(1) * b(2) - a(2) * b(1) 
281: END FUNCTION cross 
282:     end subroutine compute_py_values261:     end subroutine compute_py_values
283: 262: 
284:     subroutine compute_contact(elli, ellj, lambda_c, F)263:     subroutine compute_contact(elli, ellj, lambda_c, F)
285:         use commons, only: paramonovpbcx, paramonovpbcy, paramonovpbcz, boxlx, boxly, boxlz, &264:         use commons, only: paramonovpbcx, paramonovpbcy, paramonovpbcz, boxlx, boxly, boxlz, &
286:             & paramonovcutoff, pcutoff, pycfthresh, pycoldfusion265:             & paramonovcutoff, pcutoff, pycfthresh, pycoldfusion
287: 266: 
288:         implicit none267:         implicit none
289:         type(py_ellipsoid), intent(inout) :: elli, ellj268:         type(py_ellipsoid), intent(inout) :: elli, ellj
290:         double precision, intent(out) :: lambda_c, F269:         double precision, intent(out) :: lambda_c, F
291: 270: 
292:         integer :: image_list_offset, xyz, o1, o2, image_num, i271:         integer :: image_list_offset, xyz, o1, o2, image_num
293:         double precision :: pbc_values(3), image_rj(8, 3), rij(3), rij_p2, rijmod272:         double precision :: pbc_values(3), image_rj(8, 3), rij(3), rij_p2, rijmod
294:         double precision :: this_lambda_c, d_R, this_d_R, this_F273:         double precision :: this_lambda_c, d_R, this_d_R, this_F
295: 274: 
296:         ! make sure images are initially set properly275:         ! make sure images are initially set properly
297:         above_cutoff = .false.276:         above_cutoff = .false.
298:         elli%r(:) = elli%site%r(:)277:         elli%r(:) = elli%site%r(:)
299:         ellj%r(:) = ellj%site%r(:)278:         ellj%r(:) = ellj%site%r(:)
300:         rij(:) = elli%r(:) - ellj%r(:)279:         rij(:) = elli%r(:) - ellj%r(:)
301:         elli%mat_dot_r(:) = matmul(elli%shapemat, elli%r) ! seems that introducing this line solves the discontinuity for PBCs! 
302:         ellj%mat_dot_r(:) = matmul(ellj%shapemat, ellj%r)280:         ellj%mat_dot_r(:) = matmul(ellj%shapemat, ellj%r)
303: !                    write(*,*) 'initial distances', rij(:) 
304: 281: 
305:         if(paramonovpbcx .or. paramonovpbcy .or. paramonovpbcz) then282:         if(paramonovpbcx .or. paramonovpbcy .or. paramonovpbcz) then
306:             ! set up the pbc matrix283:             ! set up the pbc matrix
307: !        WRITE(*,*) 'in compute_contact, boxlx=', boxlx, pcutoff 
308:             pbc_values(:) = 0.d0284:             pbc_values(:) = 0.d0
309:             if(paramonovpbcx) pbc_values(1) = boxlx285:             if(paramonovpbcx) pbc_values(1) = boxlx
310:             if(paramonovpbcy) pbc_values(2) = boxly286:             if(paramonovpbcy) pbc_values(2) = boxly
311:             if(paramonovpbcz) pbc_values(3) = boxlz287:             if(paramonovpbcz) pbc_values(3) = boxlz
 288: 
312:             do xyz = 1, 3289:             do xyz = 1, 3
313:                 if(pbc_values(xyz) /= 0.d0 .and. abs(rij(xyz)) > pbc_values(xyz) / 2.d0) then290:                 if(pbc_values(xyz) /= 0.d0 .and. abs(rij(xyz)) > pbc_values(xyz) / 2.d0) &
314: !                    write(*,*) 'distance, box length',xyz, rij(xyz), pbc_values(xyz) 
315:                     ellj%r(xyz) = ellj%r(xyz) + sign(pbc_values(xyz), rij(xyz))291:                     ellj%r(xyz) = ellj%r(xyz) + sign(pbc_values(xyz), rij(xyz))
316:                 end if    
317:             end do292:             end do
318: 293: 
319:             ! redo dependent values, but don't change any orientational stuff294:             ! redo dependent values, but don't change any orientational stuff
320:             rij(:) = elli%r(:) - ellj%r(:)295:             rij(:) = elli%r(:) - ellj%r(:)
321:             ellj%mat_dot_r(:) = matmul(ellj%shapemat, ellj%r)296:             ellj%mat_dot_r(:) = matmul(ellj%shapemat, ellj%r)
322: !                    write(*,*) 'modified distances', rij(:) 
323: 297: 
324:         end if298:         end if
325: 299: 
326:         rijmod = sqrt(dot_product(rij, rij))300:         rijmod = sqrt(dot_product(rij, rij))
327: 301: 
328:         if(paramonovcutoff .and. rijmod - elli%long_semiaxis - ellj%long_semiaxis > pcutoff) then302:         if(paramonovcutoff .and. rijmod - elli%long_semiaxis - ellj%long_semiaxis > pcutoff) then
329:             ! above cutoff303:             ! above cutoff
330:             lambda_c = 0.d0304:             lambda_c = 0.d0
331:             F = 10.d0305:             F = -1.d0
332:             above_cutoff = .true.306:             above_cutoff = .true.
333:         else307:         else
334:             call polyecf(elli%shapemat_inv, ellj%shapemat_inv, elli%shapemat_inv_cofactors, &308:             call polyecf(elli%shapemat_inv, ellj%shapemat_inv, elli%shapemat_inv_cofactors, &
335:                 & ellj%shapemat_inv_cofactors, elli%shapemat_inv_det, ellj%shapemat_inv_det, & 309:                 & ellj%shapemat_inv_cofactors, elli%shapemat_inv_det, ellj%shapemat_inv_det, & 
336:                 & rij, lambda_c, F)310:                 & rij, lambda_c, F)
337:                 if(F < pycfthresh) pycoldfusion = .true.311:                 if(F < pycfthresh) pycoldfusion = .true.
338:         end if312:         end if
339:         ! if ECF falls below given value, then flag for cold fusion of the ellipsoids313:         ! if ECF falls below given value, then flag for cold fusion of the ellipsoids
340: 314: 
341:     end subroutine compute_contact315:     end subroutine compute_contact
347:         implicit none321:         implicit none
348:         double precision, parameter :: pi = 3.14159265358979323846322:         double precision, parameter :: pi = 3.14159265358979323846
349:         type(py_molecule), target, intent(inout) :: mol323:         type(py_molecule), target, intent(inout) :: mol
350:         type(py_site), pointer :: site324:         type(py_site), pointer :: site
351:         double precision, intent(inout) :: r(3), p(3)   ! position and orientation vectors325:         double precision, intent(inout) :: r(3), p(3)   ! position and orientation vectors
352:         logical, intent(in) :: gtest    ! if gradients are required326:         logical, intent(in) :: gtest    ! if gradients are required
353: 327: 
354:         integer :: i328:         integer :: i
355:         double precision :: pmod    ! modulus of p329:         double precision :: pmod    ! modulus of p
356: 330: 
357: !        write(*,*) 'before update', r(1), nint(r(1) / boxlx) 
358:  
359:         ! put molecule back in PBC box331:         ! put molecule back in PBC box
360:         if(paramonovpbcx) r(1) = r(1) - boxlx * nint(r(1) / boxlx)332:         if(paramonovpbcx) r(1) = r(1) - boxlx * nint(r(1) / boxlx)
361:         if(paramonovpbcy) r(2) = r(2) - boxly * nint(r(2) / boxly)333:         if(paramonovpbcy) r(2) = r(2) - boxly * nint(r(2) / boxly)
362:         if(paramonovpbcz) r(3) = r(3) - boxlz * nint(r(3) / boxlz)334:         if(paramonovpbcz) r(3) = r(3) - boxlz * nint(r(3) / boxlz)
363:         mol%r(:) = r(:)335:         mol%r(:) = r(:)
364: !        write(*,*) 'after update', r(1), nint(r(1) / boxlx)336: 
365:         ! make sure that 0 < |p| < 2*pi337:         ! make sure that 0 < |p| < 2*pi
366:         pmod = sqrt(dot_product(p, p))338:         pmod = sqrt(dot_product(p, p))
367:         if(pmod > 2 * pi) p(:) = p(:) / pmod * mod(pmod, 2 * pi)339:         if(pmod > 2 * pi) p(:) = p(:) / pmod * mod(pmod, 2 * pi)
368:         mol%p(:) = p(:)340:         mol%p(:) = p(:)
369: 341: 
370:         ! recompute the rotation matrices and derivatives342:         ! recompute the rotation matrices and derivatives
371:         call rmdrvt(mol%p, mol%rotmat, &343:         call rmdrvt(mol%p, mol%rotmat, &
372:             & mol%rotmat_deriv(1, :, :), mol%rotmat_deriv(2, :, :), mol%rotmat_deriv(3, :, :), gtest)344:             & mol%rotmat_deriv(1, :, :), mol%rotmat_deriv(2, :, :), mol%rotmat_deriv(3, :, :), gtest)
373: 345: 
374:         ! loop over sites346:         ! loop over sites
783:     implicit none755:     implicit none
784:     logical, intent(in) :: gtest756:     logical, intent(in) :: gtest
785:     double precision, intent(inout) :: x(3 * natoms)757:     double precision, intent(inout) :: x(3 * natoms)
786:     double precision, intent(out) :: energy, grad(3 * natoms)758:     double precision, intent(out) :: energy, grad(3 * natoms)
787: 759: 
788:     type(py_molecule), pointer :: moli, molj760:     type(py_molecule), pointer :: moli, molj
789:     type(py_site), pointer :: sitei, sitej761:     type(py_site), pointer :: sitei, sitej
790:     integer :: num_mols762:     integer :: num_mols
791:     integer :: moli_index, molj_index, sitei_index, sitej_index, i1763:     integer :: moli_index, molj_index, sitei_index, sitej_index, i1
792:     double precision :: energy_contrib, grad_contrib(12), energy_contrib2, grad_contrib2764:     double precision :: energy_contrib, grad_contrib(12), energy_contrib2, grad_contrib2
793:     double precision :: t_rep(6), t_att(6) 
794: 765: 
795:     num_mols = size(molecules)766:     num_mols = size(molecules)
796: 767: 
797:     energy = 0.d0768:     energy = 0.d0
798:     vt(:) = 0.d0769:     vt(:) = 0.d0
799:     if(gtest) grad(:) = 0.d0770:     if(gtest) grad(:) = 0.d0
800: 771: 
801:     ! update the values for all the molecules772:     ! update the values for all the molecules
802:     do moli_index = 1, num_mols773:     do moli_index = 1, num_mols
803:         moli => molecules(moli_index)774:         moli => molecules(moli_index)
829:                 sitei => moli%sites(sitei_index)800:                 sitei => moli%sites(sitei_index)
830: 801: 
831:                 ! loop over sites in inner molecule802:                 ! loop over sites in inner molecule
832:                 do sitej_index = 1, size(molj%sites)803:                 do sitej_index = 1, size(molj%sites)
833:                     sitej => molj%sites(sitej_index)804:                     sitej => molj%sites(sitej_index)
834: 805: 
835:                     energy_contrib = 0.d0806:                     energy_contrib = 0.d0
836:                     grad_contrib(:) = 0.d0807:                     grad_contrib(:) = 0.d0
837: 808: 
838:                     ! compute the energy and gradient for this pair of sites809:                     ! compute the energy and gradient for this pair of sites
839:                     call pairwise_py(sitei, sitej, grad_contrib, energy_contrib, gtest, t_rep, t_att)810:                     call pairwise_py(sitei, sitej, grad_contrib, energy_contrib, gtest)
840: 811: 
841:                     ! add the energy to total and pairwise812:                     ! add the energy to total and pairwise
842:                     energy = energy + energy_contrib813:                     energy = energy + energy_contrib
843:                     vt(moli_index) = vt(moli_index) + energy_contrib814:                     vt(moli_index) = vt(moli_index) + energy_contrib
844:                     vt(molj_index) = vt(molj_index) + energy_contrib815:                     vt(molj_index) = vt(molj_index) + energy_contrib
845: 816: 
846:                     if(gtest) then  ! add gradient contributions817:                     if(gtest) then  ! add gradient contributions
847:                         grad(moli%r_index: moli%r_index + 2) = grad(moli%r_index: moli%r_index + 2) + grad_contrib(1: 3)818:                         grad(moli%r_index: moli%r_index + 2) = grad(moli%r_index: moli%r_index + 2) + grad_contrib(1: 3)
848:                         grad(molj%r_index: molj%r_index + 2) = grad(molj%r_index: molj%r_index + 2) + grad_contrib(4: 6)819:                         grad(molj%r_index: molj%r_index + 2) = grad(molj%r_index: molj%r_index + 2) + grad_contrib(4: 6)
849:                         grad(moli%p_index: moli%p_index + 2) = grad(moli%p_index: moli%p_index + 2) + grad_contrib(7: 9)820:                         grad(moli%p_index: moli%p_index + 2) = grad(moli%p_index: moli%p_index + 2) + grad_contrib(7: 9)
970:     min_vt = minval(vt)941:     min_vt = minval(vt)
971: 942: 
972:     ! update ellipsoids before taking any steps943:     ! update ellipsoids before taking any steps
973:     do moli_index = 1, num_mols944:     do moli_index = 1, num_mols
974:         moli => molecules(moli_index)945:         moli => molecules(moli_index)
975:         r => x(moli%r_index: moli%r_index + 2)946:         r => x(moli%r_index: moli%r_index + 2)
976:         p => x(moli%p_index: moli%p_index + 2)947:         p => x(moli%p_index: moli%p_index + 2)
977: 948: 
978:         if(.not. percolatet) then949:         if(.not. percolatet) then
979:             if(vec_len(moli%r) > radius) then950:             if(vec_len(moli%r) > radius) then
980:                 write(myunit, *) 'py_takestep> initial coord outside container, bringing in'951:                 write(myunit, *) 'py_takestep> initial coord outside container, brining in'
981:                 r(:) = r(:) - sqrt(radius) * nint(r(:) / sqrt(radius))952:                 r(:) = r(:) - sqrt(radius) * nint(r(:) / sqrt(radius))
982:             end if953:             end if
983:         end if954:         end if
984: 955: 
985:         call update_py_molecule(moli, r, p, .false.)956:         call update_py_molecule(moli, r, p, .false.)
986:     end do957:     end do
987: 958: 
988:     do moli_index = 1, num_mols959:     do moli_index = 1, num_mols
989:         moli => molecules(moli_index)960:         moli => molecules(moli_index)
990:         r => x(moli%r_index: moli%r_index + 2)961:         r => x(moli%r_index: moli%r_index + 2)


r30609/rigidmd.f90 2016-07-06 15:36:11.815714865 +0100 r30608/rigidmd.f90 2016-07-06 15:36:13.595738939 +0100
  1: !  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/rigidmd.f90' in revision 30608
  2: ! svn info  
  3: ! $Date: 2011-10-16 20:53:28 +0100 (Sun, 16 Oct 2011) $ 
  4: ! $Rev: 61 $ 
  5: ! 
  6:  
  7: !     ------------------------------------------------------------------------------------ 
  8:  
  9:       MODULE MDCOMMONS 
 10:  
 11:       IMPLICIT NONE 
 12:  
 13:       SAVE 
 14:  
 15: ! N = max number of particles? used for setting array sizes 
 16: ! NN = max number of total neighbors.  used for setting array sizes 
 17: ! NMOL = number of molecules 
 18: ! NSITE = number of sites per molecule?  normally 3.  still used, but 
 19: !         sometimes assumed to be three 
 20: ! LIST = list of lists of all neighbors of each molecule. 
 21: ! POINT = holds the starting indices for each molecule for navigating 
 22: !         LIST.  e.g.  LIST( POINT(I) ) is the first neighbor of 
 23: !         molecule I 
 24: ! MASS = total mass of molecule 
 25: ! MST = holds the XYZ coordinates of the 3 beads of the molecule 
 26: !       relative to the center of mass 
 27: ! R0  = hold the xyz coords of the center of mass of the molecules at 
 28: !       last save 
 29: ! JMI = moment of inertia of the molecule 
 30: ! R = array of xyz coords of the center of mass of the molecules 
 31: ! QTRN = array of quaternions defining the angular orientation of the 
 32: !        molecules 
 33: ! P = the conjugate momentum associated with QTRN 
 34: ! V = velocity of the CoM of the molecules 
 35: !     1/MASS*(the momentum conjugate to R) 
 36: ! W = angular velocity of the molecules body-frame 
 37: ! F = forces on the molecules  
 38: ! T = torques on the molecules 
 39: ! PEP = potential energy of the molecules? 
 40: ! DELT = step size? 
 41: ! 
 42: ! BOXL = box length 
 43: ! KEP = twice? the kinetic energy of each particle ? 
 44: ! PEP = potential energy of each particle 
 45: ! SUMEPT = sum of total energy over all time steps times DELT 
 46: ! RCUT  = RCUT is the cutoff in the Lennard Jones potential. 
 47: ! RCUTSQ = RCUT*RCUT 
 48: ! RCUT2 = Maximum molecule to molecule separation at which to calculate 
 49: !         forces.  WARNING, this is measured from 
 50: !         CoM to CoM, so it will exclude sites which are roughly RCUT2-1.33 appart. 
 51: ! RCUT2SQ = RCUT2*RCUT2 
 52: ! RLIST = Maximum molecule to molecule separation at which to build up neighbor 
 53: !         list 
 54: ! RLSTSQ = RLIST**2 
 55: ! UPDATE = if true then update neighbor list 
 56: ! PHI = an array keeping track of the angular deviation of each molecule 
 57: !       from it's initial position (at time t=NEQ).  It is calulated as 
 58: !       a sum over time of W*DELT, so it is slightly approximate. 
 59: ! 
 60: ! S = S is the additional degree of freedom in the Nose-Poincare thermostat 
 61: ! PS = the conjugate momentum associated with S 
 62: ! HNOT = the value of the hamiltonitan at time t=0.  Could also have been called 
 63: !        HNAUGHT or H0   
 64:       INTEGER, PARAMETER :: N = 100, NN = N * 250 
 65:       INTEGER          :: NMOL, NSITE=2, NTST, LIST(NN), POINT(N), G 
 66:       DOUBLE PRECISION :: R(N,3), QTRN(N,4), V(N,3), W(N,3), F(N,3), T(N,3), T4(N,4), P(N,4), PHI(N,3)  
 67:       DOUBLE PRECISION :: MST(3,3), R0(N,3)  
 68:       DOUBLE PRECISION :: KEP(N), PEP(N), SUMEPT(N) 
 69:       DOUBLE PRECISION :: BOXL, RCUT, RCUTSQ, RCUT2, RCUT2SQ, EPS4, CNSTA, CNSTB, DELT, HALFDT, MASS, JMI(3) 
 70:       DOUBLE PRECISION :: RLIST, RLSTSQ, TMPFIX 
 71:       LOGICAL          :: UPDATE 
 72:       !these are used only in the NVT RUN 
 73:       DOUBLE PRECISION :: FCTMQT=1, S=1, PS=1, HNOT=1 
 74:  
 75:       END MODULE MDCOMMONS 
 76:  
 77: !     ---------------------------------------------------------------------------------------------- 
 78:  
 79:       SUBROUTINE FORQUE(PE, VRL) 
 80:  
 81: !     CALCULATES FORCES AND TORQUES on the molecules FROM SITE-SITE POTENTIAL  
 82:  
 83:       USE MDCOMMONS 
 84:       use py_module 
 85:       use commons, only : pyt, boxlx, boxly, boxlz 
 86:       IMPLICIT NONE 
 87:   
 88:       INTEGER          :: I, J, INDX, K, JBEG, JEND, JNEB, NLIST  
 89:       DOUBLE PRECISION :: RSITE(NTST,3), FSITE(NTST,3) 
 90:       !DOUBLE PRECISION :: Q(4), Q2Q3, Q1Q4, Q2Q4, Q1Q3, Q3Q4, Q1Q2, RMX(3,3) 
 91:       !DOUBLE PRECISION :: DR(3), DSS(3), DCMCM, DSS2, R2, R6, R12!, SHIFT(3) 
 92:       DOUBLE PRECISION :: DR(3), DCMCM  
 93:       DOUBLE PRECISION :: PE, VRL, AA(3), PERIODIC_OFFSET(NMOL,3) 
 94:       double precision :: energy_contrib, grad_contrib(12),INDX1,INDX2,t_rep(6),t_att(6) 
 95:       type(py_molecule), pointer :: moli, molj 
 96:       type(py_site), pointer :: sitei, sitej 
 97:       integer :: moli_index, molj_index, sitei_index, sitej_index, j1, j2, FRAMENUM 
 98:  
 99: !     INITIALIZE 
100:  
101:       PE         = 0.D0 ! total potential energy 
102:       VRL        = 0.D0 ! accumulates force*position. used to calulcate pressure 
103:       PEP(:)     = 0.D0 ! potential energy of the molecules 
104:       F(:,:)     = 0.D0 ! force on the molecules in space-fixed frame 
105:       FSITE(:,:) = 0.D0 ! force on each site in space-fixed frame 
106:       T(:,:)     = 0.D0 ! torque on the molecules 
107:       FRAMENUM = 1 
108:       PERIODIC_OFFSET(:,:) = 0.0D0       
109: !     OBTAIN THE SITE POSITIONS IN THE SPACE-FIXED FRAME FROM THE 
110: !     CENTRE-OF-MASS POSITION & THE SITE POSITIONS IN THE BODY-FIXED 
111: !     FRAME USING THE ROTATION MATRIX 
112:  
113: !      CALL GET_RSITE( RSITE ) 
114:  
115: ! get orientation from quaternions to angle-axis 
116: !        WRITE(*,*) 'before apply_periodic' 
117: !        WRITE(*,*) r(:,:) 
118:  
119: !      CALL APPLY_PERIODIC( PERIODIC_OFFSET, 1 ) 
120: !        WRITE(*,*) 'after apply_periodic' 
121: !        WRITE(*,*) r(:,:) 
122:  
123:     ! update the values for all the molecules 
124:     do moli_index = 1, nmol 
125:         CALL QTRN2AA( QTRN(moli_index,:), AA ) 
126:         moli => molecules(moli_index) 
127:         call update_py_molecule(moli, r(moli_index,:), & 
128:             & AA, .true.) 
129: !        WRITE(*,*) moli_index 
130:     end do 
131: IF (UPDATE) THEN  !if update neighbor list 
132:  
133: !     SAVE CURRENT POSITIONS 
134:  
135:          CALL SVCNFG() 
136:  
137: !     CONSTRUCT NEIGHBOUR LIST AND CALCULATE FORCES AND TORQUES 
138:  
139:          NLIST = 0 
140:  
141: !         DO I = 1, NMOL - 1 
142:  
143: !            POINT(I) = NLIST + 1 
144:  
145: !            DO J =  I + 1, NMOL 
146:  
147:   do moli_index = 1, nmol - 1 
148:         moli => molecules(moli_index) 
149:         I = moli_index 
150:              POINT(I) = NLIST + 1 
151:         
152:         ! inner loop over molecules 
153:         do molj_index = moli_index + 1, nmol 
154:             molj => molecules(molj_index) 
155:             J = molj_index 
156:                ! CoM separation  
157:                DR(:)  = R(I,:) - R(J,:) - NINT((R(I,:) - R(J,:))/BOXLX)*BOXLX 
158:                DCMCM  = DOT_PRODUCT(DR,DR) 
159:  
160: !               IF (DCMCM < RLSTSQ ) THEN 
161:  
162:                   NLIST       = NLIST + 1 
163:                   LIST(NLIST) = J 
164:  
165: !       REMOVE THIS CHECK IF MAXNAB IS APPROPRIATE 
166:  
167:                   IF (NLIST == NN) STOP 'list too small' 
168:  
169:                 IF(PYT) THEN 
170:                     ! loop over sites in outer molecule 
171:                   do sitei_index = 1, size(moli%sites) 
172:                       sitei => moli%sites(sitei_index) 
173:  
174:                       ! loop over sites in inner molecule 
175:                       do sitej_index = 1, size(molj%sites) 
176:                           sitej => molj%sites(sitej_index) 
177:  
178:                           energy_contrib = 0.d0 
179:                           grad_contrib(:) = 0.d0 
180:  
181:                           ! compute the energy and gradient for this pair of sites 
182:                           call pairwise_py(sitei, sitej, grad_contrib, energy_contrib, .true., t_rep, t_att) 
183:                           ! add the energy to total  
184:                           PE = PE + energy_contrib 
185:                           PEP(I) = PEP(I) + energy_contrib 
186:                           PEP(J) = PEP(J) + energy_contrib 
187:                            DO K = 1, 3 ! accumulate forces and torques 
188:                                 F(I,k) = F(I,k) - grad_contrib(k) 
189:                                 F(J,k) = F(J,k) - grad_contrib(k+3) 
190:                                 T(I,k) = T(I,k) - t_rep(k) - t_att(k)  
191:                                 T(J,k) = T(J,k) - t_rep(3+k) - t_att(3+k) 
192:                            END DO 
193:  
194:                       end do 
195:                   end do 
196:  
197:                 ELSE ! LWOTP model 
198:                  CALL MOLMOL_FORQUE(I, J, RSITE, FSITE, DR, DCMCM, PE, VRL) 
199:                 END IF 
200: !               ENDIF 
201:  
202:             ENDDO 
203:                     
204:          ENDDO  
205:  
206:         POINT(NMOL) = NLIST + 1  
207: !        WRITE(*,*) NMOL, POINT(NMOL), NLIST 
208: ELSE ! if .not. UPDATE use current neighbor list to find forces 
209:  
210:   do moli_index = 1, nmol - 1 
211:         moli => molecules(moli_index) 
212:         I = moli_index 
213: !         DO I = 1, NMOL - 1 
214:  
215: !       USE THE LIST TO FIND THE NEIGHBOURS     
216:  
217:             JBEG = POINT(I) 
218:             JEND = POINT(I+1) - 1 
219:  
220: !       CHECK THAT MOLECULE I HAS NEIGHBOURS 
221:  
222:       IF (JBEG <= JEND) THEN 
223: !                WRITE(*,*) JBEG, JEND, NMOL, BOXL 
224:         do JNEB = JBEG, JEND 
225: !               DO JNEB = JBEG, JEND 
226:  
227:                   J = LIST(JNEB) 
228:  
229:                   DR(:)  = R(I,:) - R(J,:) - NINT((R(I,:) - R(J,:))/BOXL)*BOXL 
230:                   DCMCM  = DOT_PRODUCT(DR,DR) 
231:                 IF(PYT) THEN 
232:                         molj_index = J 
233:                         molj => molecules(molj_index) 
234:                             ! loop over sites in outer molecule 
235:                   do sitei_index = 1, size(moli%sites) 
236:                               sitei => moli%sites(sitei_index) 
237:  
238:                       ! loop over sites in inner molecule 
239:                       do sitej_index = 1, size(molj%sites) 
240:                           sitej => molj%sites(sitej_index) 
241:  
242:                           energy_contrib = 0.d0 
243:                           grad_contrib(:) = 0.d0 
244:  
245:                           ! compute the energy and gradient for this pair of sites 
246:                           call pairwise_py(sitei, sitej, grad_contrib, energy_contrib, .true., t_rep, t_att) 
247:                           ! add the energy to total  
248:                           PE = PE + energy_contrib 
249:                           PEP(I) = PEP(I) + energy_contrib 
250:                           PEP(J) = PEP(J) + energy_contrib 
251:  
252:                            DO K = 1, 3 ! accumulate forces and torques 
253:                                 F(I,k) = F(I,k) - grad_contrib(k) 
254:                                 F(J,k) = F(J,k) - grad_contrib(k+3) 
255: !                                T(I,k) = T(I,k) - grad_contrib(6+k) 
256: !                                T(J,k) = T(J,k) - grad_contrib(9+k) 
257:                                 T(I,k) = T(I,k) - t_rep(k) - t_att(k)  
258:                                 T(J,k) = T(J,k) - t_rep(3+k) - t_att(3+k) 
259:                            END DO 
260:                       end do 
261:                   end do 
262:  
263:  
264:                 ELSE ! LWOTP model 
265:  
266:                   CALL MOLMOL_FORQUE(I, J, RSITE, FSITE, DR, DCMCM, PE, VRL) 
267:                 END IF 
268:                ENDDO 
269:  
270:             ENDIF 
271:  
272:          ENDDO 
273:  
274:       ENDIF 
275:  
276: !     MULTIPLY RESULTS BY ENERGY FACTORS  
277:       IF(.NOT.PYT) THEN 
278:               FSITE = FSITE * EPS4 * 2.D0 
279:               PE    = PE * EPS4 
280:               VRL   = VRL * EPS4 * 2.D0 
281:       END IF 
282: !     CONVERT SITE FORCES TO MOLECULAR FORCE AND TORQUE 
283:  
284:       DO I = 1, NMOL 
285:          DO J = 1, NSITE 
286:  
287:             INDX = (I - 1) * NSITE + J 
288:  
289:           IF(PYT) THEN 
290: !              CALL AA2QTRN(T(I,:), 
291: !            T(I,1) = T(I,1) - grad_contrib(8) + grad_contrib(9) 
292: !            T(I,2) = T(I,2) - grad_contrib(9) + grad_contrib(7) 
293: !            T(I,3) = T(I,3) - grad_contrib(7) + grad_contrib(8)  
294:           ELSE 
295:             DO K = 1, 3 
296:                F(I,K) = F(I,K) + FSITE(INDX,K) 
297:             ENDDO 
298:             T(I,1) = T(I,1) + RSITE(INDX,2)*FSITE(INDX,3) - RSITE(INDX,3)*FSITE(INDX,2) 
299:             T(I,2) = T(I,2) + RSITE(INDX,3)*FSITE(INDX,1) - RSITE(INDX,1)*FSITE(INDX,3) 
300:             T(I,3) = T(I,3) + RSITE(INDX,1)*FSITE(INDX,2) - RSITE(INDX,2)*FSITE(INDX,1)  
301:           END IF 
302:          ENDDO 
303:        
304:       ENDDO 
305:  
306:       END SUBROUTINE FORQUE 
307:   
308: !     ---------------------------------------------------------------------------------------------- 
309:  
310:       SUBROUTINE DEFMOL() 
311:  
312: !     DEFINE THE MOLECULAR GEOMETRY WITH THE CENTRE-OF-MASS AT THE  
313: !     ORIGIN & CALCULATE THE TOTAL MASS & THE ORINCIPAL OF MOMENTS OF  
314: !     INERTIA 
315:  
316:       USE MDCOMMONS, ONLY : MST, JMI, MASS, MST, RCUT2, RCUT2SQ, RCUT, NSITE 
317:  
318:       IMPLICIT NONE 
319:  
320:       INTEGER :: I 
321:       DOUBLE PRECISION :: MS(NSITE), PI, DIST(NSITE), DELRC 
322:  
323:       PI       = 4.D0 * DATAN(1.D0) 
324:  
325:       MST(1,1) = 0.D0 
326:       MST(1,2) = - 2.D0 * DSIN(7.D0*PI/24.D0) / 3.D0 
327:       MST(1,3) = 0.D0 
328:  
329:       MST(2,1) = DCOS(7.D0*PI/24.D0) 
330:       MST(2,2) = DSIN(7.D0*PI/24.D0) / 3.D0 
331:       MST(2,3) = 0.D0 
332:  
333:       MST(3,1) = - DCOS(7.D0*PI/24.D0) 
334:       MST(3,2) = DSIN(7.D0*PI/24.D0) / 3.D0 
335:       MST(3,3) = 0.D0 
336:  
337:       DO I = 1, NSITE 
338:  
339:          DIST(I) = DSQRT(DOT_PRODUCT(MST(I,:),MST(I,:))) 
340:  
341:       ENDDO 
342:  
343:       DELRC    = 2.D0 * MAXVAL(DIST)  
344:       RCUT2    = RCUT + DELRC 
345:       RCUT2SQ  = RCUT2 * RCUT2 
346:   
347:       MS(1)    = 1.D0/3.D0 
348:       MS(2)    = 1.D0/3.D0 
349:       MS(3)    = 1.D0/3.D0 
350:  
351:       MASS     = MS(1) + MS(2) + MS(3) 
352:  
353:       JMI(1)  = MS(1)*MST(1,2)**2+MS(2)*MST(2,2)**2+MS(3)*MST(3,2)**2 
354:       JMI(2)  = MS(1)*MST(1,1)**2+MS(2)*MST(2,1)**2+MS(3)*MST(3,1)**2 
355:       JMI(3)  = MS(1)*MST(1,1)**2+MS(2)*MST(2,1)**2+MS(3)*MST(3,1)**2   & 
356:                  + MS(1)*MST(1,2)**2+MS(2)*MST(2,2)**2+MS(3)*MST(3,2)**2 
357:       JMI(:) = 1.0D0 
358:       MST(:,:) = 1.0D0 
359: !        WRITE(*,*) 'JMI=', JMI(:), MASS, MST, NSITE 
360: !STOP 
361:       END SUBROUTINE DEFMOL    
362:         
363: !     ---------------------------------------------------------------------------------------------- 
364:  
365:       SUBROUTINE EFM (EFMT, TIME) 
366:  
367:       ! calculate EFMT, the energy fluctuation metric: see  
368:       ! de Sauza and Wales J. Chem. Phys. 123 134504 2005  
369:       ! http://dx.doi.org/10.1063/1.2035080 
370:       ! calculate energy averaged over time and molecules. 
371:       ! calculate variance over molecules of time averaged energy 
372:  
373:       USE MDCOMMONS, ONLY : SUMEPT, NMOL, DELT, PEP, KEP 
374:  
375:       IMPLICIT NONE 
376:  
377:       INTEGER          :: I 
378:       DOUBLE PRECISION :: AVET, EFMT, EP, TIME, DBLE 
379:  
380:       TIME = TIME + DELT 
381:       AVET = 0.D0 ! energy averaged over time and molecules <<E>_T>_N 
382:       EFMT = 0.D0 ! variance of average energy over molecules? <(<E>_T-AVET)^2>_N 
383:  
384:       DO I = 1, NMOL 
385:  
386:          EP        = 0.5D0 * KEP(I) + PEP(I) 
387:          SUMEPT(I) = SUMEPT(I) + EP * DELT 
388:          AVET      = AVET + SUMEPT(I)  
389:  
390:       ENDDO 
391:  
392:       AVET = AVET / (TIME * DBLE(NMOL)) 
393:  
394:       DO I = 1, NMOL 
395:  
396:          EFMT = EFMT + (SUMEPT(I) / TIME - AVET) ** 2.D0 
397:  
398:       ENDDO 
399:  
400:       EFMT = EFMT / DBLE(NMOL) 
401:  
402:       END SUBROUTINE EFM 
403:  
404: !     ---------------------------------------------------------------------------------------------- 
405:  
406:       SUBROUTINE RIGIDMD_SC () 
407:  
408: ! initialize the CoM positions on a simple cubic lattice 
409:  
410:       USE MDCOMMONS, only : R, NMOL 
411:       USE COMMONS, only : boxlx 
412:  
413:       IMPLICIT NONE 
414:  
415:       INTEGER          :: I, NUC, IX, IY, IZ, M!, K 
416:       DOUBLE PRECISION :: BOXLH, UCL, UCLH, C(3), DBLE 
417:  
418: !     NUMBER OF UNIT CELLS 
419:       BOXLH = 0.5D0 * BOXLX 
420:  
421:       I = 1 
422:  
423:       DO WHILE (I ** 3 < NMOL) 
424:          I = I + 1 
425:       ENDDO 
426:  
427:       NUC   = I 
428:  
429: !     UNIT CELL LENGTH 
430:  
431:       UCL     = BOXLX / DBLE(NUC) 
432:       UCLH    = 0.5D0 * UCL 
433:  
434: !     CONSTRUCT THE LATTICE FROM THE UNIT CELL 
435:  
436:       M = 0 
437:  
438:       DO IZ = 1, NUC 
439:  
440:          C(3) = IZ * UCL - UCLH - BOXLH 
441:  
442:          DO IY = 1, NUC 
443:  
444:             C(2) = IY * UCL - UCLH - BOXLH 
445:  
446:             DO IX = 1, NUC 
447:  
448:                C(1) = IX * UCL - UCLH - BOXLH 
449:  
450:                M = M + 1 
451:                R(M,:) = C(:) 
452:  
453:             ENDDO 
454:  
455:          ENDDO 
456:  
457:       ENDDO 
458:  
459:       END SUBROUTINE RIGIDMD_SC 
460:  
461: !     ---------------------------------------------------------------------------------------------- 
462:  
463:       SUBROUTINE EQCON() 
464:  
465: ! load positions, orientations, linear velocities, and angular 
466: ! velocities from files 
467:  
468:       USE MDCOMMONS , only : R, QTRN, V, W, NMOL 
469:  
470:       IMPLICIT NONE 
471:  
472:       INTEGER :: I 
473:  
474:       OPEN (UNIT=13, FILE='initpos1.dat', STATUS='UNKNOWN') 
475:       OPEN (UNIT=14, FILE='initortn1.dat', STATUS='UNKNOWN') 
476:       OPEN (UNIT=15, FILE='initlinvel1.dat', STATUS='UNKNOWN') 
477:       OPEN (UNIT=16, FILE='initangvel1.dat', STATUS='UNKNOWN') 
478:  
479:       DO I = 1, NMOL   
480:          
481:          READ(13, *) R(I,1), R(I,2), R(I,3) 
482:          !READ(14, *) QTRN(I,3), QTRN(I,2), QTRN(I,4), QTRN(I,1) 
483:          READ(14, *) QTRN(I,1), QTRN(I,2), QTRN(I,3), QTRN(I,4) 
484:          READ(15, *) V(I,1), V(I,2), V(I,3) 
485:          READ(16, *) W(I,1), W(I,2), W(I,3) 
486:          !QTRN(I,3) = -QTRN(I,3) 
487:  
488:       ENDDO 
489:  
490:        
491:       CLOSE (13) 
492:       CLOSE (14) 
493:       CLOSE (15) 
494:       CLOSE (16) 
495:  
496:       END SUBROUTINE EQCON 
497:  
498: !     ---------------------------------------------------------------------------------------------- 
499:  
500:       SUBROUTINE LOAD_POS_ORTN() 
501:  
502: ! load positions, orientations, linear velocities, and angular 
503: ! velocities from files 
504:  
505:       USE MDCOMMONS , only : R, QTRN, NMOL 
506:  
507:       IMPLICIT NONE 
508:  
509:       INTEGER :: I 
510:  
511:       OPEN (UNIT=13, FILE='initpos1.dat', STATUS='UNKNOWN') 
512:       OPEN (UNIT=14, FILE='initortn1.dat', STATUS='UNKNOWN') 
513:       !OPEN (UNIT=15, FILE='initlinvel1.dat', STATUS='UNKNOWN') 
514:       !OPEN (UNIT=16, FILE='initangvel1.dat', STATUS='UNKNOWN') 
515:         WRITE(*,*) 'NMOL=', NMOL 
516:       DO I = 1, NMOL   
517:          
518:          READ(13, *) R(I,1), R(I,2), R(I,3) 
519:          !READ(14, *) QTRN(I,3), QTRN(I,2), QTRN(I,4), QTRN(I,1) 
520:          READ(14, *) QTRN(I,1), QTRN(I,2), QTRN(I,3), QTRN(I,4) 
521:          !READ(15, *) V(I,1), V(I,2), V(I,3) 
522:          !READ(16, *) W(I,1), W(I,2), W(I,3) 
523:          !QTRN(I,3) = -QTRN(I,3) 
524:          WRITE(*,*) "Read molecule ", I 
525:          WRITE(*,*) R(I,1), R(I,2), R(I,3) 
526:       ENDDO 
527:  
528:        
529:       CLOSE (13) 
530:       CLOSE (14) 
531:       !CLOSE (15) 
532:       !CLOSE (16) 
533:  
534:       END SUBROUTINE LOAD_POS_ORTN 
535:  
536: !     ---------------------------------------------------------------------------------------------- 
537:  
538:       SUBROUTINE LOAD_LIN_ANG_VEL() 
539:  
540: ! load positions, orientations, linear velocities, and angular 
541: ! velocities from files 
542:  
543:       USE MDCOMMONS , only : V, W, NMOL 
544:  
545:       IMPLICIT NONE 
546:  
547:       INTEGER :: I 
548:  
549:       !OPEN (UNIT=13, FILE='initpos1.dat', STATUS='UNKNOWN') 
550:       !OPEN (UNIT=14, FILE='initortn1.dat', STATUS='UNKNOWN') 
551:       OPEN (UNIT=15, FILE='initlinvel1.dat', STATUS='UNKNOWN') 
552:       OPEN (UNIT=16, FILE='initangvel1.dat', STATUS='UNKNOWN') 
553:  
554:       DO I = 1, NMOL   
555:          
556:          !READ(13, *) R(I,1), R(I,2), R(I,3) 
557:          !READ(14, *) QTRN(I,3), QTRN(I,2), QTRN(I,4), QTRN(I,1) 
558:          !READ(14, *) QTRN(I,1), QTRN(I,2), QTRN(I,3), QTRN(I,4) 
559:          READ(15, *) V(I,1), V(I,2), V(I,3) 
560:          READ(16, *) W(I,1), W(I,2), W(I,3) 
561:          !QTRN(I,3) = -QTRN(I,3) 
562:  
563:       ENDDO 
564:  
565:        
566:       CLOSE (13) 
567:       CLOSE (14) 
568:       !CLOSE (15) 
569:       !CLOSE (16) 
570:  
571:       END SUBROUTINE LOAD_LIN_ANG_VEL 
572:  
573: !     ---------------------------------------------------------------------------------------------- 
574:  
575:       SUBROUTINE INTVEL(TMPINT) 
576:  
577: ! initialize linear velocities and angular velocity 
578:  
579:       USE MDCOMMONS, only : MASS, JMI, V, NMOL, W 
580:  
581:       IMPLICIT NONE 
582:  
583:       INTEGER          :: I 
584:       DOUBLE PRECISION :: TMPINT, VMAG, WMAG, SUMV(3), AVV(3), E(3), DBLE 
585:        
586:       VMAG    = DSQRT(3.D0 * TMPINT / MASS) 
587:       WMAG = DSQRT(3.D0 * TMPINT / (JMI(1) + JMI(2) + JMI(3))) 
588:       SUMV(:) = 0.D0 
589:  
590:       DO I = 1, NMOL 
591:  
592:          CALL RNDVCT(E) 
593:  
594:          V(I,:) = VMAG * E(:) 
595:          SUMV(:) = SUMV(:) + V(I,:) 
596:             
597:       ENDDO  
598:  
599:       AVV(:) = SUMV(:) / DBLE(NMOL) 
600:       
601:       DO I = 1, NMOL 
602:  
603:          V(I,:) = V(I,:) - AVV(:) 
604:  
605:          CALL RNDVCT(E) 
606:          W(I,:) = WMAG * E(:) 
607:  
608:       ENDDO 
609:  
610:       END SUBROUTINE INTVEL  
611:  
612: !     ---------------------------------------------------------------------------------------------- 
613:  
614:       SUBROUTINE INTORN() 
615:  
616: ! initialize orientations 
617:  
618:       USE MDCOMMONS, only : NMOL, QTRN 
619:  
620:       IMPLICIT NONE 
621:  
622:       INTEGER          :: I 
623:       DOUBLE PRECISION :: NORM 
624:  
625: !     SET QUATERNIONS WITH RANDOM ORIENTATIONS 
626:        
627:       DO I =1, NMOL 
628:  
629:          QTRN(I,1) = 0.10 
630:          QTRN(I,2) = -0.415D0 
631:          QTRN(I,3) = -0.246D0 
632:          QTRN(I,4) = -0.976D0 
633:  
634:          NORM   = DSQRT(DOT_PRODUCT(QTRN(I,:),QTRN(I,:))) 
635:  
636:          QTRN(I,:) = QTRN(I,:) / NORM 
637:   
638:       ENDDO 
639:  
640:       END SUBROUTINE INTORN 
641:  
642: !     ---------------------------------------------------------------------------------------------- 
643:  
644:       SUBROUTINE RNDVCT(E) 
645:  
646: !     CHOOSE A RANDOM UNIT VECTOR IN SPACE 
647:  
648:       IMPLICIT NONE 
649:  
650:       DOUBLE PRECISION :: E(3), DUMMY, XI, XISQ, XI1, XI2 
651:  
652:       XISQ = 2.D0 
653:  
654:       DO WHILE (XISQ > 1.D0) 
655:  
656:          CALL RANDOM_NUMBER (DUMMY) 
657:  
658:          XI1  = DUMMY * 2.D0 - 1.D0 
659:  
660:          CALL RANDOM_NUMBER (DUMMY) 
661:  
662:          XI2  = DUMMY * 2.D0 - 1.D0 
663:          XISQ = XI1 * XI1 + XI2 * XI2 
664:  
665:       ENDDO 
666:  
667:       XI    = DSQRT(1.D0 - XISQ) 
668:  
669:       E(1) = 2.D0 * XI1 * XI 
670:       E(2) = 2.D0 * XI2 * XI 
671:       E(3) = 1.D0 - 2.D0 * XISQ 
672:  
673:       END SUBROUTINE RNDVCT 
674:  
675: !     ---------------------------------------------------------------------------------------------- 
676:  
677:       SUBROUTINE CHECK () 
678:  
679:       USE MDCOMMONS, ONLY : NMOL, R, R0, RCUT2, RLIST, UPDATE 
680:  
681: !     DECIDES WHETHER THE LIST NEEDS TO BE RECONSTRUCTED FROM THE COORDINATES AT LAST UPDATE 
682:  
683: !     LOGICAL  UPDATE    IF TRUE THE LIST IS UPDATED 
684:  
685: !     CHECK IS CALLED TO SET UPDATE BEFORE EVERY CALL TO FORQUE 
686:  
687:       IMPLICIT NONE  
688:  
689:       INTEGER          :: I, K 
690:       DOUBLE PRECISION :: DISPMX 
691:  
692: !     CALCULATE MAXIMUM DISPLACEMENT SINCE LAST UPDATE 
693:  
694:       DISPMX = 0.D0 
695:  
696:       DO I = 1, NMOL 
697:  
698:          DO K = 1, 3 
699:  
700:             DISPMX = MAX(ABS(R(I,K) - R0(I,K)), DISPMX) 
701:  
702:          ENDDO 
703: !        WRITE(*,*) 'DISPMX=', DISPMX 
704:       ENDDO 
705:  
706: !     A CONSERVATIVE TEST OF THE LIST SKIN CROSSING 
707:  
708:       DISPMX = 2.D0 * DSQRT(3.D0 * DISPMX ** 2) 
709:  
710:       UPDATE = (DISPMX > (RLIST - RCUT2)) 
711: !      WRITE(*,*) 'UPDATE', UPDATE 
712:       END SUBROUTINE CHECK  
713:  
714: !     ---------------------------------------------------------------------------------------------- 
715:  
716:       SUBROUTINE SVCNFG() 
717:  
718: ! save config 
719:  
720:       USE MDCOMMONS, only : R, R0, NMOL 
721:         
722: !     SAVE IS CALLED WHENEVER THE NEW VERLET LIST IS CONSTRUCTED 
723:  
724:       IMPLICIT NONE 
725:  
726:       INTEGER :: I, K 
727:  
728:       DO I = 1, NMOL 
729:  
730:          DO K = 1, 3 
731:  
732:            R0(I,K) = R(I,K) 
733:  
734:          ENDDO 
735:  
736:       ENDDO 
737:  
738:       END SUBROUTINE SVCNFG 
739:  
740: !     ---------------------------------------------------------------------------------------------- 
741:  
742:       SUBROUTINE RUN_INFO() 
743:         
744: !     print some info about the run.  time elapsed, etc. 
745:  
746:       !USE MDCOMMONS 
747:  
748:       IMPLICIT NONE 
749:  
750:       real :: timearray(2), timesecs 
751:  
752:       OPEN (UNIT = 40, FILE = 'runinfo.dat', STATUS = 'UNKNOWN') 
753:  
754: !      call etime( timearray, timesecs ) 
755:  
756:       write(40,*) "time_elapsed ", timesecs, timearray 
757:  
758:       close(40) 
759:  
760:       END SUBROUTINE RUN_INFO 
761:  
762: !     ---------------------------------------------------------------------------------------------- 
763:  
764:       SUBROUTINE GET_RSITE( RSITE ) 
765:         
766: !     GET the positions of the sites with respect to the CoM 
767: !     i.e. RSITE_fixed_space(I,:) = RSITE(I,:) + R(I,:) 
768: !     use MST and the orientations QTRN 
769:  
770:       USE MDCOMMONS, only : Nmol, QTRN, NSITE, MST, NTST 
771:  
772:       IMPLICIT NONE 
773:  
774:       DOUBLE PRECISION :: RSITE(NTST,3) 
775:       DOUBLE PRECISION :: Q(4), RMX(3,3) 
776:       integer :: INDX, I, J 
777:  
778:       DO I = 1, NMOL 
779:  
780: !     CONSTRUCT THE ROTATION MATRIX RMX FROM THE QUATERNION 
781:  
782:          Q    = QTRN(I,:) 
783:  
784:          CALL CONSTRUCT_RMX( RMX, Q) 
785:  
786:          DO J = 1, NSITE 
787:  
788:             INDX = (I - 1) * NSITE + J 
789:  
790:             !the site coordinate in fixed space frame 
791:             RSITE(INDX,:) = MATMUL(RMX,MST(J,:)) 
792:      
793:          ENDDO 
794:  
795:       ENDDO 
796:  
797:       END SUBROUTINE GET_RSITE 
798:  
799: !     ---------------------------------------------------------------------------------------------- 
800:  
801:       SUBROUTINE MOLMOL_FORQUE(I, J, RSITE, FSITE, DR, DCMCM, PE, VRL) 
802:         
803: !     GET the positions of the sites with respect to the CoM 
804: !     i.e. RSITE_fixed_space(I,:) = RSITE(I,:) + R(I,:) 
805: !     use MST and the orientations QTRN 
806:  
807:       USE MDCOMMONS, only : RCUT2SQ, NSITE, RCUTSQ, PEP, NTST, CNSTA, CNSTB 
808:  
809:       IMPLICIT NONE 
810:  
811:       INTEGER          :: I, J, J1, J2, INDX1, INDX2, K 
812:       DOUBLE PRECISION :: RSITE(NTST,3), FSITE(NTST,3) 
813:       !DOUBLE PRECISION :: Q(4), Q2Q3, Q1Q4, Q2Q4, Q1Q3, Q3Q4, Q1Q2, RMX(3,3) 
814:       DOUBLE PRECISION :: DR(3), DSS(3), DCMCM, DSS2, R2, R6, R12!, SHIFT(3) 
815:       DOUBLE PRECISION :: PE, VJ1J2, FJ1J2, VRL 
816:  
817:       IF (DCMCM < RCUT2SQ) THEN 
818:  
819:         INDX1 = (I - 1) * NSITE 
820:         INDX2 = (J - 1) * NSITE 
821:  
822:         DO J1 = 1, NSITE 
823:  
824:           DO J2 = 1, NSITE 
825:  
826:             ! DSS = distance between the sites 
827:             DSS(:) = DR(:) + RSITE(INDX1+J1,:) - RSITE(INDX2+J2,:) 
828:             DSS2   = DOT_PRODUCT(DSS,DSS) 
829:  
830:             if (DSS2 < RCUTSQ ) THEN 
831:               R2     = 1.D0 / DSS2 
832:               R6     = R2 * R2 * R2 
833:               R12    = R6 ** 2 
834:               VJ1J2  = R12 - R6 + CNSTA * DSS2 + CNSTB 
835:               PE     = PE + VJ1J2 
836:               PEP(I) = PEP(I) + 2.D0 * VJ1J2 
837:               PEP(J) = PEP(J) + 2.D0 * VJ1J2 
838:               FJ1J2  = (6.D0 * R12 - 3.D0 * R6)/DSS2 - CNSTA !half the force 
839:               VRL    = VRL + FJ1J2 * DSS2 !used to calculate the pressure 
840:  
841:               DO K = 1, 3 
842:  
843:                 FSITE(INDX1+J1,K) = FSITE(INDX1+J1,K) + FJ1J2*DSS(K) 
844:                 FSITE(INDX2+J2,K) = FSITE(INDX2+J2,K) - FJ1J2*DSS(K) 
845:  
846:               ENDDO 
847:  
848:             ENDIF 
849:  
850:           ENDDO 
851:  
852:         ENDDO 
853:  
854:       ENDIF 
855:  
856:       END SUBROUTINE MOLMOL_FORQUE 
857:  
858: !     ---------------------------------------------------------------------------------------------- 
859:  
860:       SUBROUTINE CONSTRUCT_MX( MX, Q ) 
861:  
862: !     CONSTRUCT THE MATRIX MX FROM THE QUATERNION Q 
863: ! MX = the matrix representation of Q 
864: !  Q1 -Q2 -Q3 -Q4  
865: !  Q2  Q1 -Q4  Q3 
866: !  Q3  Q4  Q1 -Q2 
867: !  Q4 -Q3  Q2  Q1 
868: !        OR, with indices from 0 to 3 
869: !  Q0 -Q1 -Q2 -Q3  
870: !  Q1  Q0 -Q3  Q2 
871: !  Q2  Q3  Q0 -Q1 
872: !  Q3 -Q2  Q1  Q0 
873:  
874:  
875:       !USE MDCOMMONS 
876:  
877:       IMPLICIT NONE 
878:  
879:       DOUBLE PRECISION, INTENT(OUT) :: MX(4,4) 
880:       DOUBLE PRECISION, INTENT(IN) :: Q(4) 
881:       INTEGER :: K 
882:  
883:       DO K = 1, 4 
884:         MX(K,K) = Q(1) 
885:       ENDDO 
886:  
887:       MX(1,2:4) = -Q(2:4) 
888:       MX(2:4,1) = Q(2:4) 
889:       MX(2,3)   = -Q(4) 
890:       MX(2,4)   = Q(3) 
891:       MX(3,2)   = Q(4) 
892:       MX(3,4)   = -Q(2) 
893:       MX(4,2)   = -Q(3) 
894:       MX(4,3)   = Q(2) 
895:  
896:       END SUBROUTINE CONSTRUCT_MX 
897:  
898: !     ---------------------------------------------------------------------------------------------- 
899:  
900:       SUBROUTINE CONSTRUCT_RMX( RMX, Q ) 
901:  
902:       !construct the rotation matrix which rotates a 3-vector from the body 
903:       !fixed frame to the space fixed frame 
904:       !uses the fact that the quaternion is normalized 
905:       !1=a^2+b^2+c^2+d^2 
906:  
907:       !USE MDCOMMONS 
908:  
909:       IMPLICIT NONE 
910:  
911:       DOUBLE PRECISION, INTENT(OUT) :: RMX(3,3) 
912:       DOUBLE PRECISION, INTENT(IN) :: Q(4) 
913:       DOUBLE PRECISION :: Q2Q3, Q1Q4, Q2Q4, Q1Q3, Q3Q4, Q1Q2  
914:  
915:       Q2Q3 = Q(2)*Q(3) 
916:       Q1Q4 = Q(1)*Q(4) 
917:       Q2Q4 = Q(2)*Q(4) 
918:       Q1Q3 = Q(1)*Q(3) 
919:       Q3Q4 = Q(3)*Q(4) 
920:       Q1Q2 = Q(1)*Q(2) 
921:  
922:       RMX(1,1) = 0.5D0 - Q(3)*Q(3) - Q(4)*Q(4) 
923:       RMX(2,2) = 0.5D0 - Q(2)*Q(2) - Q(4)*Q(4) 
924:       RMX(3,3) = 0.5D0 - Q(2)*Q(2) - Q(3)*Q(3) 
925:       RMX(1,2) = Q2Q3 - Q1Q4 
926:       RMX(2,1) = Q2Q3 + Q1Q4 
927:       RMX(1,3) = Q2Q4 + Q1Q3 
928:       RMX(3,1) = Q2Q4 - Q1Q3 
929:       RMX(2,3) = Q3Q4 - Q1Q2 
930:       RMX(3,2) = Q3Q4 + Q1Q2 
931:  
932:       RMX      = 2.D0*RMX 
933:  
934:       END SUBROUTINE CONSTRUCT_RMX 
935:  
936: !     ---------------------------------------------------------------------------------------------- 
937:  
938:       SUBROUTINE CONSTRUCT_RMXT( RMXT, Q ) 
939:  
940:       !construct the rotation matrix which rotates a 3-vector from the space 
941:       !fixed frame to the body frame 
942:       !uses the fact that the quaternion is normalized 
943:       !1=a^2+b^2+c^2+d^2 
944:  
945:       !USE MDCOMMONS 
946:  
947:       IMPLICIT NONE 
948:  
949:       DOUBLE PRECISION, INTENT(OUT) :: RMXT(3,3) 
950:       DOUBLE PRECISION, INTENT(IN) :: Q(4) 
951:       DOUBLE PRECISION :: Q2Q3, Q1Q4, Q2Q4, Q1Q3, Q3Q4, Q1Q2  
952:  
953:       Q2Q3 = Q(2)*Q(3) 
954:       Q1Q4 = Q(1)*Q(4) 
955:       Q2Q4 = Q(2)*Q(4) 
956:       Q1Q3 = Q(1)*Q(3) 
957:       Q3Q4 = Q(3)*Q(4) 
958:       Q1Q2 = Q(1)*Q(2) 
959:  
960:       RMXT(1,1) = 0.5D0 - Q(3)*Q(3) - Q(4)*Q(4) 
961:       RMXT(2,2) = 0.5D0 - Q(2)*Q(2) - Q(4)*Q(4) 
962:       RMXT(3,3) = 0.5D0 - Q(2)*Q(2) - Q(3)*Q(3) 
963:       RMXT(1,2) = Q2Q3 + Q1Q4 
964:       RMXT(2,1) = Q2Q3 - Q1Q4 
965:       RMXT(1,3) = Q2Q4 - Q1Q3 
966:       RMXT(3,1) = Q2Q4 + Q1Q3 
967:       RMXT(2,3) = Q3Q4 + Q1Q2 
968:       RMXT(3,2) = Q3Q4 - Q1Q2 
969:  
970:       RMXT      = 2.D0*RMXT 
971:  
972:       END SUBROUTINE CONSTRUCT_RMXT 
973:  
974: !     ---------------------------------------------------------------------------------------------- 
975:  
976:       SUBROUTINE QTRN2AA( Q, V ) 
977:  
978:       !convert from quaternion to angle axis representation of a 
979:       !rotation 
980:  
981:       IMPLICIT NONE 
982:  
983:       DOUBLE PRECISION, INTENT(OUT) :: V(3) 
984:       DOUBLE PRECISION, INTENT(IN) :: Q(4) 
985:       DOUBLE PRECISION :: THETA, D 
986:  
987:       THETA = 2.D0*ACOS(Q(1)) 
988:       D     = DSQRT(1.D0-Q(1)*Q(1)) 
989:       V(:)  = THETA * Q(2:4)/D 
990:  
991:       END SUBROUTINE QTRN2AA 
992:  
993: !     ---------------------------------------------------------------------------------------------- 
994:  
995:       SUBROUTINE AA2QTRN( Q, V ) 
996:  
997:       !convert from angle axis to quaternion representation of a rotation 
998:  
999:       IMPLICIT NONE 
1000:  
1001:       DOUBLE PRECISION, INTENT(IN) :: V(3) 
1002:       DOUBLE PRECISION, INTENT(OUT) :: Q(4) 
1003:       DOUBLE PRECISION :: THETA 
1004:  
1005:       THETA = SQRT( V(1)* V(1)+ V(2)* V(2)+ V(3)* V(3) ) 
1006:       Q(1)    = COS( THETA / 2.D0 ) 
1007:       Q(2:4) = SQRT( 1-Q(1)*Q(1) ) * V(:) / THETA 
1008:  
1009:       END SUBROUTINE AA2QTRN 
1010:  
1011:  
1012: !     ---------------------------------------------------------------------------------------------- 
1013:  
1014:       SUBROUTINE PfromW( P, W, JMI, Q, S ) 
1015:  
1016:       !convert from angle axis to quaternion representation of a rotation 
1017:  
1018:       IMPLICIT NONE 
1019:  
1020:       DOUBLE PRECISION, INTENT(IN) :: W(3), JMI(3), Q(4), S 
1021:       DOUBLE PRECISION, INTENT(OUT) :: P(4) 
1022:       DOUBLE PRECISION :: W4(4), MX(4,4) 
1023:  
1024:       ! intitialize P, the angular momentum.  P is virtual, W is real 
1025:       CALL CONSTRUCT_MX( MX, Q ) 
1026:       W4(1)  = 0.D0 
1027:       W4(2)  = 2.D0*W(1)*JMI(1) 
1028:       W4(3)  = 2.D0*W(2)*JMI(2) 
1029:       W4(4)  = 2.D0*W(3)*JMI(3) 
1030:       P(:) = MATMUL(MX*S,W4) 
1031:  
1032:  
1033:       END SUBROUTINE PfromW 
1034:  
1035: !     ---------------------------------------------------------------------------------------------- 
1036:  
1037:       SUBROUTINE WfromP( P, W, JMI, Q, S ) 
1038:  
1039:       !convert from angle axis to quaternion representation of a rotation 
1040:  
1041:       IMPLICIT NONE 
1042:  
1043:       DOUBLE PRECISION, INTENT(IN) :: P(4), JMI(3), Q(4), S 
1044:       DOUBLE PRECISION, INTENT(OUT) :: W(3) 
1045:       DOUBLE PRECISION :: W4(4), MX(4,4) 
1046:  
1047:       CALL CONSTRUCT_MX( MX, Q ) 
1048:       W4     = MATMUL(TRANSPOSE(MX),(P(:)/S)) 
1049:       W(1) = 0.5D0*W4(2)/JMI(1) 
1050:       W(2) = 0.5D0*W4(3)/JMI(2) 
1051:       W(3) = 0.5D0*W4(4)/JMI(3) 
1052:  
1053:       END SUBROUTINE WfromP 
1054:  
1055: ! 
1056: ! svn info  
1057: ! $Date: 2011-10-16 20:53:28 +0100 (Sun, 16 Oct 2011) $ 
1058: ! $Rev: 61 $ 
1059: ! 
1060:  
1061:   
1062: !     ---------------------------------------------------------------------------------------------- 
1063:   
1064:       SUBROUTINE MOVE (ISTEP, NEQ, PE, VRL, TKE, RKE) 
1065:  
1066: !     propagator: H. Okumura, S. G. Itoh, and Y. Okamoto, J. Chem. Phys. 126, 084103 (2007). 
1067:       !note: move operates on whole molecules, except through the subroutine 
1068:       !FORQUE 
1069:  
1070:       USE MDCOMMONS 
1071:       USE COMMONS, only : PYT 
1072:  
1073:       IMPLICIT NONE 
1074:  
1075:       ! PHIT = This parameter is called ZETA in the above reference.  As far as 
1076:       !        I can tell it is not anything physical.  one can define it 
1077:       !        component by component like  
1078:       !        PHIT(J) = DOT_PRODUCT( P(I,:) * MX(:,J) ) / (4 * JMI(J) * S) 
1079:       !        where MX is the matrix representation of QTRN(I,:) as defined 
1080:       !        below 
1081:       INTEGER          :: ISTEP, NEQ, I 
1082:       !DOUBLE PRECISION :: MX(4,4), W4(4), RMXT(3,3)  
1083:       DOUBLE PRECISION :: RMXT(3,3)  
1084:       DOUBLE PRECISION :: PE, VRL, TKE, RKE, RKET 
1085:       DOUBLE PRECISION :: Q(4),T4dummy(4) 
1086:  
1087:       RKE = 0.D0 
1088:  
1089:       ! S = S is the additional degree of freedom in the Nose-Poincare thermostat 
1090:       ! MQ = the artificial mass associated with S 
1091:       ! PS = the conjugate momentum associated with S 
1092:  
1093:  
1094:  
1095:       ! STEP 1: update S, PS 
1096:       CALL NVESTEP1( S, PS, FCTMQT ) 
1097:  
1098:       ! STEP 2: update PS, V from F, and P from T4 
1099:       CALL NVESTEP2(PE) 
1100:  
1101:       ! STEP 3: update QTRN, P, PS with contributions from principle axis 3 
1102:       CALL NVESTEP3() 
1103:        
1104:       ! STEP 4: update QTRN, P, PS with contributions from principle axis 2 
1105:       CALL NVESTEP4() 
1106:  
1107:       ! STEP 5: update QTRN, P, PS with contributions from principle axis 1.  Also update R from V 
1108:       CALL NVESTEP5() 
1109:  
1110:       ! STEP 6: update QTRN, P, PS with contributions from principle axis 2 
1111:       ! STEP 6 is exactly the same as STEP 4 
1112:       CALL NVESTEP4() 
1113:  
1114:       ! STEP 7: update QTRN, P, PS with contributions from principle axis 3 
1115:       ! STEP 7 is exactly the same as STEP 3 
1116:       CALL NVESTEP3() 
1117:  
1118:        
1119:       !question: why is FORQUE called here????? 
1120:       CALL CHECK () 
1121:       CALL FORQUE (PE, VRL) 
1122:       ! update T4 
1123:       DO I = 1, NMOL 
1124:          Q = QTRN(I,:) 
1125:          CALL CONSTRUCT_RMXT( RMXT, Q ) 
1126:          T4(I,1)   = 0.D0 
1127: !       IF(PYT) THEN 
1128:          
1129: !         CALL AA2QTRN(T4dummy(:),T(I,:)) 
1130: !         T4(I,2:4) = MATMUL(RMXT,T4dummy(2:4)) 
1131: !         T4(I,1)   = 0.D0 
1132: !       ELSE 
1133:          T4(I,2:4) = MATMUL(RMXT,T(I,:)) 
1134: !       END IF 
1135:       ENDDO 
1136:  
1137:       ! STEP 8: update PS, V from F, and P from T4 
1138:       ! STEP 8 is exactly the same as STEP 2 
1139:       CALL NVESTEP2(PE) 
1140:  
1141:       ! STEP 9: update S, PS 
1142:       ! STEP 9 is exactly the same as STEP 1 
1143:       !RKET, RKE, PHI, etc.???? 
1144:       CALL NVESTEP1( S, PS, FCTMQT ) 
1145:  
1146:  
1147:  
1148:       ! calculate W from P (uses MX) 
1149:       ! update KEP, TKE, RKET, RKE, PHI 
1150:       TKE    = 0.D0 
1151:       DO I = 1, NMOL 
1152:          KEP(I) = MASS * (V(I,1)*V(I,1) + V(I,2)*V(I,2) + V(I,3)*V(I,3))/(S*S) !here V is virtual velocity 
1153:           
1154:          TKE = TKE + KEP(I)  
1155:  
1156: !        ADVANCE ANGULAR VELOCITY 
1157:  
1158:          Q    = QTRN(I,:) 
1159:  
1160:          !NOTE: if I incorporate the W update into STEP 8 I can avoid 
1161:          !this extra call to CONSTRUCT_MX 
1162:          !CALL CONSTRUCT_MX( MX, Q ) 
1163:          !W4     = MATMUL(TRANSPOSE(MX),(P(I,:)/S)) 
1164:          !W(I,1) = 0.5D0*W4(2)/JMI(1) 
1165:          !W(I,2) = 0.5D0*W4(3)/JMI(2) 
1166:          !W(I,3) = 0.5D0*W4(4)/JMI(3) 
1167:          CALL WfromP( P(I,:), W(I,:), JMI, Q, S ) 
1168:  
1169:          RKET   = JMI(1)*W(I,1)*W(I,1) + JMI(2)*W(I,2)*W(I,2) + JMI(3)*W(I,3)*W(I,3) 
1170:          KEP(I) = KEP(I) + RKET 
1171:          RKE    = RKE + RKET 
1172: !         RKE = 0 
1173:          IF (ISTEP > NEQ) THEN 
1174:  
1175:             PHI(I,:) = PHI(I,:) + W(I,:) * DELT 
1176:  
1177:          ENDIF 
1178:  
1179:       ENDDO  
1180:  
1181:       TKE = 0.5D0 * TKE 
1182:       RKE = 0.5D0 * RKE 
1183:  
1184:  
1185:  
1186:       END SUBROUTINE MOVE 
1187: !     ---------------------------------------------------------------------------------------------- 
1188:   
1189:       SUBROUTINE NVTMOVE (ISTEP, NEQ, PE, VRL, TKE, RKE) 
1190:  
1191: !     propagator: H. Okumura, S. G. Itoh, and Y. Okamoto, J. Chem. Phys. 126, 084103 (2007). 
1192:       !note: move operates on whole molecules, except through the subroutine 
1193:       !FORQUE 
1194:  
1195:       USE MDCOMMONS 
1196:       USE COMMONS, only : PYT 
1197:  
1198:       IMPLICIT NONE 
1199:  
1200:       ! PHIT = This parameter is called ZETA in the above reference.  As far as 
1201:       !        I can tell it is not anything physical.  one can define it 
1202:       !        component by component like  
1203:       !        PHIT(J) = DOT_PRODUCT( P(I,:) * MX(:,J) ) / (4 * JMI(J) * S) 
1204:       !        where MX is the matrix representation of QTRN(I,:) as defined 
1205:       !        below 
1206:       INTEGER          :: ISTEP, NEQ, I 
1207:       !DOUBLE PRECISION :: MX(4,4), W4(4), RMXT(3,3)  
1208:       DOUBLE PRECISION :: RMXT(3,3)  
1209:       DOUBLE PRECISION :: PE, VRL, TKE, RKE, RKET 
1210:       DOUBLE PRECISION :: Q(4),T4dummy(4) 
1211:  
1212:       RKE = 0.D0 
1213:  
1214:       ! S = S is the additional degree of freedom in the Nose-Poincare thermostat 
1215:       ! MQ = the artificial mass associated with S 
1216:       ! PS = the conjugate momentum associated with S 
1217:  
1218:  
1219:  
1220:       ! STEP 1: update S, PS 
1221:       CALL NVTSTEP1( S, PS, FCTMQT ) 
1222:  
1223:       ! STEP 2: update PS, V from F, and P from T4 
1224:       CALL NVTSTEP2(PE) 
1225:  
1226:       ! STEP 3: update QTRN, P, PS with contributions from principle axis 3 
1227:       CALL NVTSTEP3() 
1228:        
1229:       ! STEP 4: update QTRN, P, PS with contributions from principle axis 2 
1230:       CALL NVTSTEP4() 
1231:  
1232:       ! STEP 5: update QTRN, P, PS with contributions from principle axis 1.  Also update R from V 
1233:       CALL NVTSTEP5() 
1234:  
1235:       ! STEP 6: update QTRN, P, PS with contributions from principle axis 2 
1236:       ! STEP 6 is exactly the same as STEP 4 
1237:       CALL NVTSTEP4() 
1238:  
1239:       ! STEP 7: update QTRN, P, PS with contributions from principle axis 3 
1240:       ! STEP 7 is exactly the same as STEP 3 
1241:       CALL NVTSTEP3() 
1242:  
1243:        
1244:       !question: why is FORQUE called here????? 
1245:       CALL CHECK () 
1246:       CALL FORQUE (PE, VRL) 
1247:       ! update T4 
1248:       DO I = 1, NMOL 
1249:          Q = QTRN(I,:) 
1250:          CALL CONSTRUCT_RMXT( RMXT, Q ) 
1251:          T4(I,1)   = 0.D0 
1252: !       IF(PYT) THEN 
1253:          
1254: !         CALL AA2QTRN(T4dummy(:),T(I,:)) 
1255: !         T4(I,2:4) = MATMUL(RMXT,T4dummy(2:4)) 
1256: !         T4(I,1)   = 0.D0 
1257: !       ELSE 
1258:          T4(I,2:4) = MATMUL(RMXT,T(I,:)) 
1259: !       END IF 
1260:       ENDDO 
1261:  
1262:       ! STEP 8: update PS, V from F, and P from T4 
1263:       ! STEP 8 is exactly the same as STEP 2 
1264:       CALL NVTSTEP2(PE) 
1265:  
1266:       ! STEP 9: update S, PS 
1267:       ! STEP 9 is exactly the same as STEP 1 
1268:       !RKET, RKE, PHI, etc.???? 
1269:       CALL NVTSTEP1( S, PS, FCTMQT ) 
1270:  
1271:  
1272:  
1273:       ! calculate W from P (uses MX) 
1274:       ! update KEP, TKE, RKET, RKE, PHI 
1275:       TKE    = 0.D0 
1276:       DO I = 1, NMOL 
1277:          KEP(I) = MASS * (V(I,1)*V(I,1) + V(I,2)*V(I,2) + V(I,3)*V(I,3))/(S*S) !here V is virtual velocity 
1278:           
1279:          TKE = TKE + KEP(I)  
1280:  
1281: !        ADVANCE ANGULAR VELOCITY 
1282:  
1283:          Q    = QTRN(I,:) 
1284:  
1285:          !NOTE: if I incorporate the W update into STEP 8 I can avoid 
1286:          !this extra call to CONSTRUCT_MX 
1287:          !CALL CONSTRUCT_MX( MX, Q ) 
1288:          !W4     = MATMUL(TRANSPOSE(MX),(P(I,:)/S)) 
1289:          !W(I,1) = 0.5D0*W4(2)/JMI(1) 
1290:          !W(I,2) = 0.5D0*W4(3)/JMI(2) 
1291:          !W(I,3) = 0.5D0*W4(4)/JMI(3) 
1292:          CALL WfromP( P(I,:), W(I,:), JMI, Q, S ) 
1293:  
1294:          RKET   = JMI(1)*W(I,1)*W(I,1) + JMI(2)*W(I,2)*W(I,2) + JMI(3)*W(I,3)*W(I,3) 
1295:          KEP(I) = KEP(I) + RKET 
1296:          RKE    = RKE + RKET 
1297: !         RKE = 0 
1298:          IF (ISTEP > NEQ) THEN 
1299:  
1300:             PHI(I,:) = PHI(I,:) + W(I,:) * DELT 
1301:  
1302:          ENDIF 
1303:  
1304:       ENDDO  
1305:  
1306:       TKE = 0.5D0 * TKE 
1307:       RKE = 0.5D0 * RKE 
1308:  
1309:  
1310:  
1311:       END SUBROUTINE NVTMOVE 
1312:  
1313:       SUBROUTINE DEFMOL1() 
1314:  
1315: !     DEFINE THE MOLECULAR GEOMETRY WITH THE CENTRE-OF-MASS AT THE  
1316: !     ORIGIN & CALCULATE THE TOTAL MASS & THE ORINCIPAL OF MOMENTS OF  
1317: !     INERTIA 
1318:  
1319:       USE MDCOMMONS 
1320:  
1321:       IMPLICIT NONE 
1322:  
1323:       DOUBLE PRECISION :: MS(NSITE), PI 
1324:  
1325:       PI       = 4.D0 * DATAN(1.D0) 
1326:  
1327:       MST(1,1) = 0.D0 
1328:       MST(1,2) = - 2.D0 * DSIN(7.D0*PI/24.D0) / 3.D0 
1329:       MST(1,3) = 0.D0 
1330:  
1331:       MST(2,1) = DCOS(7.D0*PI/24.D0) 
1332:       MST(2,2) = DSIN(7.D0*PI/24.D0) / 3.D0 
1333:       MST(2,3) = 0.D0 
1334:  
1335:       MST(3,1) = - DCOS(7.D0*PI/24.D0) 
1336:       MST(3,2) = DSIN(7.D0*PI/24.D0) / 3.D0 
1337:       MST(3,3) = 0.D0 
1338:   
1339:       MS(1)    = 1.D0/3.D0 
1340:       MS(2)    = 1.D0/3.D0 
1341:       MS(3)    = 1.D0/3.D0 
1342:  
1343:       MASS     = MS(1) + MS(2) + MS(3) 
1344:  
1345:       JMI(1)  = MS(1)*MST(1,2)**2+MS(2)*MST(2,2)**2+MS(3)*MST(3,2)**2 
1346:       JMI(2)  = MS(1)*MST(1,1)**2+MS(2)*MST(2,1)**2+MS(3)*MST(3,1)**2 
1347:       JMI(3)  = MS(1)*MST(1,1)**2+MS(2)*MST(2,1)**2+MS(3)*MST(3,1)**2   & 
1348:                  + MS(1)*MST(1,2)**2+MS(2)*MST(2,2)**2+MS(3)*MST(3,2)**2 
1349:       END SUBROUTINE DEFMOL1    
1350:  
1351: !     ---------------------------------------------------------------------------------------------- 
1352:  
1353:       SUBROUTINE GET_FRAME_R_ORNT(IOstatus, AABOOL) 
1354:  
1355: ! load positions, orientations, linear velocities, and angular 
1356: ! velocities from files 
1357:  
1358:       USE MDCOMMONS, ONLY : R, QTRN, NMOL  
1359:  
1360:       IMPLICIT NONE 
1361:  
1362:       INTEGER :: I,IOstatus 
1363:       LOGICAL, INTENT(IN) :: AABOOL 
1364:       DOUBLE PRECISION :: V(3) 
1365:  
1366:       if ( AABOOL ) then 
1367:  
1368:       DO I = 1, NMOL   
1369:          
1370:          READ(13, *, IOSTAT=IOstatus) R(I,1), R(I,2), R(I,3) 
1371:          IF (IOstatus .NE. 0) return 
1372:       ENDDO 
1373:       DO I = 1, NMOL   
1374:          READ(13, *, IOSTAT=IOstatus) V(1), V(2), V(3)  
1375:          IF (IOstatus .NE. 0) return 
1376:          CALL AA2QTRN( QTRN(I,:), V ) 
1377:       ENDDO 
1378:  
1379:       ELSE 
1380:       DO I = 1, NMOL   
1381:          
1382:          READ(13, *, IOSTAT=IOstatus) R(I,1), R(I,2), R(I,3) 
1383:          IF (IOstatus .NE. 0) return 
1384:          READ(14, *, IOSTAT=IOstatus) QTRN(I,1), QTRN(I,2), QTRN(I,3), QTRN(I,4) 
1385:          IF (IOstatus .NE. 0) return 
1386:          !READ(15, *) V(I,1), V(I,2), V(I,3) 
1387:          !READ(16, *) W(I,1), W(I,2), W(I,3) 
1388:          !QTRN(I,3) = -QTRN(I,3) 
1389:  
1390:       ENDDO 
1391:  
1392:     ENDIF 
1393:  
1394:        
1395:  
1396:       END SUBROUTINE GET_FRAME_R_ORNT 
1397:  
1398: !     ---------------------------------------------------------------------------------------------- 
1399:  
1400:       SUBROUTINE SITECOORDS(RSITE) 
1401:  
1402: ! load positions, orientations, linear velocities, and angular 
1403: ! velocities from files 
1404:  
1405:       USE MDCOMMONS 
1406:  
1407:       IMPLICIT NONE 
1408:  
1409:       DOUBLE PRECISION :: Q(4), RMX(3,3) 
1410:       INTEGER          :: I, J, INDX!, INDX1, INDX2, K, JBEG, JEND, JNEB, NLIST  
1411:       DOUBLE PRECISION, INTENT(OUT) :: RSITE(NTST,3) 
1412:  
1413: !     OBTAIN THE SITE POSITIONS IN THE SPACE-FIXED FRAME FROM THE 
1414: !     CENTRE-OF-MASS POSITION & THE SITE POSITIONS IN THE BODY-FIXED 
1415: !     FRAME USING THE ROTATION MATRIX 
1416:  
1417:       DO I = 1, NMOL 
1418:  
1419: !     CONSTRUCT THE ROTATION MATRIX RMX FROM THE QUATERNION 
1420:  
1421:          Q    = QTRN(I,:) 
1422:          CALL CONSTRUCT_RMX(RMX,Q) 
1423:  
1424:          DO J = 1, NSITE 
1425:  
1426:             INDX = (I - 1) * NSITE + J 
1427:  
1428:             !the site coordinate in CoM position 
1429:             RSITE(INDX,:) = MATMUL(RMX,MST(J,:)) + R(I,:) 
1430:      
1431:          ENDDO 
1432:  
1433:       ENDDO 
1434:  
1435:       END SUBROUTINE SITECOORDS 
1436:  
1437: !     ---------------------------------------------------------------------------------------------- 
1438:  
1439:       SUBROUTINE POUT(RSITE) 
1440:  
1441:       USE MDCOMMONS 
1442:  
1443:       IMPLICIT NONE 
1444:  
1445:       DOUBLE PRECISION, INTENT(OUT) :: RSITE(NTST,3) 
1446:       INTEGER:: I 
1447:       CHARACTER(LEN=3):: NAM(3) 
1448:       NAM(1) = "N  " 
1449:       NAM(2) = "C  " 
1450:       NAM(3) = "N  " 
1451:       
1452:  
1453:       WRITE (1, *) NTST 
1454:       WRITE (1, *)  
1455:       DO I = 1,NTST 
1456:         !WRITE (1, *) "C ", RSITE(I,1), RSITE(I,2), RSITE(I,3) 
1457:         !WRITE (1, *, ADVANCE='NO' ) NAM(MODULO(I,3)+1), RSITE(I,1), RSITE(I,2), RSITE(I,3) 
1458:         WRITE (1, 900 ) NAM(MODULO(I,3)+1), RSITE(I,1), RSITE(I,2), RSITE(I,3) 
1459:         900 FORMAT (1x,A3,1x,ES20.12,1x,ES20.12,1x,ES20.12) 
1460:       ENDDO 
1461:  
1462:       END SUBROUTINE POUT 
1463:  
1464: !     ---------------------------------------------------------------------------------------------- 
1465:  
1466:       SUBROUTINE APPLY_PERIODIC( PERIODIC_OFFSET, FRAMENUM ) 
1467:  
1468:       ! apply periodic boundary conditions, but do it the same way each 
1469:       ! frame to make the visual representation more appealing. 
1470:  
1471:       USE MDCOMMONS 
1472:       USE commons, only : boxlx 
1473:  
1474:       IMPLICIT NONE 
1475:  
1476:       INTEGER :: I,J, FRAMENUM 
1477:       DOUBLE PRECISION :: PERIODIC_OFFSET(N,3) 
1478:        
1479:       IF ( FRAMENUM .EQ. 1 ) THEN 
1480:         DO I=1,NMOL 
1481:  
1482:           !apply periodic boundary conditions 
1483:           !the coordinate might be more than one box length outside of 
1484:           !the box 
1485:           DO J=1,3 
1486:             DO WHILE ( R(I,J) .GE. 0.5D0*BOXLX )  
1487:               R(I,J)=R(I,J)-BOXLX 
1488:               PERIODIC_OFFSET(I,J) = PERIODIC_OFFSET(I,J) - BOXLX 
1489:             END DO 
1490:             DO WHILE ( R(I,J) .LT. -0.5D0*BOXLX )  
1491:               R(I,J)=R(I,J)+BOXLX 
1492:               PERIODIC_OFFSET(I,J) = PERIODIC_OFFSET(I,J) + BOXLX 
1493:             END DO 
1494:           END DO 
1495:         END DO 
1496:       ELSE 
1497:         DO I=1,NMOL 
1498:           DO J=1,3 
1499:             R(I,J)=R(I,J)+ PERIODIC_OFFSET(I,J); 
1500:           END DO 
1501:         END DO 
1502:       ENDIF 
1503:  
1504:  
1505:       END SUBROUTINE APPLY_PERIODIC 
1506: ! 
1507: ! svn info  
1508: ! $Date: 2011-10-12 12:06:51 +0100 (Wed, 12 Oct 2011) $ 
1509: ! $Rev: 54 $ 
1510: ! 
1511:  
1512: ! the lewis wahnstrom model has  
1513:       ! sigma = 0.483 nm 
1514:       ! epsilon = 5.276 kJ/mol 
1515:   
1516: !     ---------------------------------------------------------------------------------------------- 
1517:  
1518: !      PROGRAM LWOTP_MD_NVE 
1519:       SUBROUTINE RIGIDMD_NVE 
1520:  
1521:       USE MDCOMMONS 
1522:       USE COMMONS, ONLY : PYT 
1523:  
1524:       IMPLICIT NONE 
1525:  
1526: ! IDUMP1: dump files.  always:  
1527: !                     enrg1.dat: KEPP, PEPP, EPP: kinetic, potential, energy per molecule 
1528: !                     tmpprs1.dat: temp, pressure? 
1529: !                     avenrg1.dat: average KE, PE, E 
1530: !                     avtmpprs1.dat: average temp, pressure? 
1531: !                     avtmptr1.dat: average tranlational, rotational temperature 
1532: ! IDUMP2: dump files.  ISTEP > NEQ 
1533: !                     pos1.dat:  xyz positions of CoM of molecules 
1534: !                     ortn1.dat: quaternions giving the orientation of molecules 
1535: !                     rot1.dat:  phi1, phi2, phi3 ???? 
1536: !                     efm1.dat:  energy fluctuation metric: time from NEQ, and variance over molecules of time averaged energy  
1537: ! IDUMP3: CALL SCLVEL.  ISTEP <= NBATH: scale linear velocities ??? 
1538: ! IDUMP4: dump files.  ISTEP > NEQ:  
1539: !                     blockav1.dat: average PE, temp, pressure 
1540: !                     lvel1.dat:    linear velocities 
1541: !                     avel1.dat:    angular velocities 
1542: ! TKE: translational kinetic energy 
1543: ! RKE: rotational kinetic energy 
1544: ! MX:  is a representation of the quaternion as a matrix in such a way 
1545: !      that quaternion addition and multiplication correspond to matrix 
1546: !      addition and matrix multiplication 
1547: !      In this way the conjugate of the quaternion corresponds to the 
1548: !      transpose of the matrix.  The fourth power of the norm of a 
1549: !      quaternion is the determinant of the corresponding matrix. 
1550: !      Complex numbers are block diagonal matrices with two 2x2 blocks. 
1551: ! EFMT: Energy Fluctuation Metric: the variance over molecules of the 
1552: !       time averaged energy 
1553: ! VRL: accumulates the forces dotted into the positions Fij.dot(Ri - Rj) which 
1554: !      is used to calculate the pressure 
1555:       INTEGER          :: NSTEP, NBATH, NEQ, ISTEP, ICNT, IDUMP1, IDUMP2, IDUMP3, IDUMP4, I!, K!, STATUS 
1556:       DOUBLE PRECISION :: RCUT6, RCUT12, TMP, TMPTR, TMPRT, RHO, VLM 
1557:       DOUBLE PRECISION :: Q(4), RMXT(3,3), MX(4,4), W4(4), T4dummy(4) 
1558:       DOUBLE PRECISION :: PRS, SUMPRS, AVPRS, PE, KE, TKE, RKE, PEPP, KEPP, EPP, VRL 
1559:       DOUBLE PRECISION :: SUME, SUMKE, SUMPE, SUMTMT, SUMTMR, SUMTMP 
1560:       DOUBLE PRECISION :: AVE, AVKE, AVPE, AVTMPT, AVTMPR, AVTMP, EFMT, TIME 
1561:       DOUBLE PRECISION :: BSUMPE, BSUMTMP, BSUMPRS, AVPEB, AVTMPB, AVPRSB, DBLE 
1562:       DOUBLE PRECISION :: MQ=0.1D0, DH=1.D0 !only NVT 
1563:       DOUBLE PRECISION :: TMPINT !initial temperature 
1564:       logical :: bool1, bool2 
1565:       !double precision :: rand 
1566:       !CALL RANDOM_SEED() 
1567:       !CALL RANDOM_number( rand ) 
1568:       !write(*,*) rand 
1569:  
1570:  
1571:       OPEN (UNIT = 1, FILE = 'parameter1.inp', STATUS = 'UNKNOWN') 
1572:  
1573: !     READ INPUT PARAMETERS 
1574:  
1575:       READ (1, *) 
1576:       READ (1, *) NMOL, NSITE 
1577:       READ (1, *)  
1578:       READ (1, *) RCUT, RLIST 
1579:       READ (1, *) 
1580:       READ (1, *) TMPFIX 
1581:       READ (1, *) 
1582:       READ (1, *) RHO 
1583:       READ (1, *) 
1584:       READ (1, *) DELT 
1585:       READ (1, *) 
1586:       READ (1, *) NSTEP, NBATH, NEQ 
1587:       READ (1, *)  
1588:       READ (1, *) IDUMP1, IDUMP2, IDUMP3, IDUMP4 
1589:   
1590:       CLOSE (1) 
1591:       TMPINT = TMPFIX !initial temperature 
1592:  
1593: !     CALCULATE FROM INPUT PARAMETERS 
1594:  
1595:       G      = 6 * NMOL - 3 
1596:       NTST   = NMOL * NSITE             ! TOTAL NUMBER OF SITES 
1597:       RCUTSQ = RCUT * RCUT 
1598:       RLSTSQ = RLIST * RLIST 
1599:       EPS4   = 4.D0  ! 4*epsilon 
1600:       RCUT6  = (1.D0 / RCUT) ** 6 
1601:       RCUT12 = RCUT6 * RCUT6 
1602:       CNSTA  = (6.D0 * RCUT12 - 3.D0 * RCUT6) / RCUTSQ 
1603:       CNSTB  = 4.D0 * RCUT6 - 7.D0 * RCUT12 
1604:       HALFDT = 0.5D0 * DELT 
1605:       FCTMQT = 0.5D0*HALFDT/MQ !used only in NVT 
1606:   
1607: !     ALLOCATE THE ARRAYS 
1608:  
1609: !      ALLOCATE (R(NMOL,3), QTRN(NMOL,4), V(NMOL,3), W(NMOL,3), F(NMOL,3), T(NMOL,3), P(NMOL,4), & 
1610: !      T4(NMOL,4), PHI(NMOL,3), FSITE(NTST,3), MST(NSITE,3), KEP(NMOL), PEP(NMOL), SUMEPT(NMOL), & 
1611: !      STAT = STATUS) 
1612:  
1613: !      IF (STATUS /= 0) STOP 'ALLOCATION PROBLEM' 
1614:   
1615: !     BOX LENGTH 
1616:  
1617:       VLM    = DBLE(NMOL) / RHO 
1618: !      BOXL   = VLM ** (1.D0/3.D0) 
1619:       BOXL = 10.0D0 
1620:       CALL DEFMOL() 
1621:  
1622:       IF ( BOXL .LE. 2.D0 + 2.D0*RCUT ) THEN 
1623:         write(*,*) "WARNING: the box is too small, molecules will be interacting with multiple images" 
1624:       ENDIF 
1625:  
1626: !      WRITE(*,*) RCUT, RCUT2, RLIST 
1627:  
1628:  
1629:       INQUIRE(FILE="initpos1.dat", EXIST=bool1 )  
1630:       INQUIRE(FILE="initortn1.dat", EXIST=bool2 )  
1631:       if ( bool1 .and. bool2 ) then 
1632:         !start from an existing configuration 
1633:         !CALL EQCON() 
1634:         write(*,*) "reading init files: initpos1.dat and initortn1.dat" 
1635:         CALL LOAD_POS_ORTN() 
1636:       else 
1637:         !     start from scratch 
1638:         !     START WITH A SIMPLE CUBIC LATTICE CONFIGURATION 
1639:         write(*,*) "starting from a simple cubic lattice" 
1640:         CALL RIGIDMD_SC() 
1641:         !     SET INITIAL ORIENTATIONS 
1642:         CALL INTORN() 
1643:       endif 
1644:  
1645: !     SET INITIAL LINEAR and angular VELICITIES 
1646:       INQUIRE(FILE="initlinvel1.dat", EXIST=bool1 )  
1647:       INQUIRE(FILE="initangvel1.dat", EXIST=bool2 )  
1648:       if ( bool1 .and. bool2 ) then 
1649:         write(*,*) "reading init files: initlinvel1.dat and initangvel1.dat" 
1650:         CALL LOAD_LIN_ANG_VEL() 
1651:       else 
1652:         write(*,*) "assigning random velocities" 
1653:         CALL INTVEL(TMPINT) 
1654:       endif 
1655:  
1656:       !for bug testing only 
1657:       IF ( .true. ) then 
1658:         OPEN (UNIT=71, FILE='pos0.dat', STATUS='UNKNOWN') 
1659:         OPEN (UNIT=72, FILE='ortn0.dat', STATUS='UNKNOWN') 
1660:         OPEN (UNIT=73, FILE='lvel0.dat', STATUS='UNKNOWN') 
1661:         OPEN (UNIT=74, FILE='avel0.dat', STATUS='UNKNOWN') 
1662:         do I=1,NMOL 
1663:           WRITE (71, *)  R(I,1), R(I,2), R(I,3) 
1664:           WRITE (72, 900)  QTRN(I,1), QTRN(I,2), QTRN(I,3), QTRN(I,4)  
1665:           WRITE (73, *)  V(I,1)/S, V(I,2)/S, V(I,3)/S 
1666:           WRITE (74, *) W(I,1), W(I,2), W(I,3) 
1667:         enddo 
1668:         close(71) 
1669:         close(72) 
1670:         close(73) 
1671:         close(74) 
1672:       endif 
1673:  
1674:  
1675: !     START SIMULATION 
1676:  
1677:       ISTEP  = 0   
1678:       ICNT   = 0 
1679:       SUME   = 0.D0 
1680:       SUMKE  = 0.D0 
1681:       SUMPE  = 0.D0 
1682:       SUMTMT = 0.D0 
1683:       SUMTMR = 0.D0 
1684:       SUMTMP = 0.D0 
1685:       SUMPRS = 0.D0 
1686:        
1687:       PHI(:,:) = 0.D0 
1688:       KE       = 0.D0 
1689:  
1690: !     CALCULATE FORCES AND TORQUES AT TIME t=0 
1691:  
1692:       UPDATE = .TRUE. 
1693:  
1694:       CALL FORQUE(PE, VRL) 
1695:  
1696:       DO I = 1, NMOL 
1697:  
1698:          !calculate the transposed rotation matrix from quaternion 
1699:          Q    = QTRN(I,:) 
1700:          CALL CONSTRUCT_RMXT( RMXT, Q ) 
1701:  
1702: !        TORQUE AT TIME t=0 IN THE BODY FRAME USING THE TRANSPOSE OF THE ROTATION MATRIX 
1703:  
1704:          T4(I,1)   = 0.D0 
1705: !         IF(PYT) THEN 
1706: !          T(I,:) = MATMUL(RMXT,T(I,:)) 
1707: !          CALL AA2QTRN(T4(I,:),T(I,:)) 
1708: !          T4(I,2:4) = MATMUL(RMXT,T4(I,2:4)) 
1709: !          T4(I,1)   = 0.D0 
1710: !          T(I,:) = MATMUL(RMXT,T(I,:)) 
1711: !          CALL AA2QTRN(T4(I,:),T(I,:)) 
1712: !          T4(I,2:4) = MATMUL(RMXT,T(I,:)) 
1713: !          T4(I,1)   = 0.D0 
1714: !          T4(I,2:4) = T(I,:) 
1715: !         ELSE 
1716:           T4(I,2:4) = MATMUL(RMXT,T(I,:)) 
1717: !         END IF 
1718:  
1719: !     CONSTRUCT THE MATRIX MX FROM THE QUATERNION Q 
1720:  
1721:          CALL CONSTRUCT_MX( MX, Q ) 
1722:  
1723:          W4(1)  = 0.D0 
1724:          W4(2)  = 2.D0*W(I,1)*JMI(1) 
1725:          W4(3)  = 2.D0*W(I,2)*JMI(2) 
1726:          W4(4)  = 2.D0*W(I,3)*JMI(3) 
1727:  
1728:          ! intitialize the momentum conjugate to angular velocity 
1729:          P(I,:) = MATMUL(MX,W4) 
1730:   
1731:          KE = KE + 0.5D0*MASS*(V(I,1)*V(I,1) + V(I,2)*V(I,2) + V(I,3)*V(I,3)) & 
1732:                  + 0.5D0*(W(I,1)*W(I,1)/JMI(1) + W(I,2)*W(I,2)/JMI(2) + W(I,3)*W(I,3)/JMI(3)) 
1733:  
1734:       ENDDO 
1735:  
1736:       HNOT = KE + PE  
1737:  
1738:       DO WHILE (ISTEP < NSTEP)             
1739:  
1740:          !check to make sure the extra variables relevant only for NVT runs are 
1741:          !constant and = 1 
1742:          !UNIT 0 should be standard error 
1743:          IF ( ABS(S-1.D0) > 1.D-7 ) write(0, *) "S != 1: ", S 
1744:          IF ( ABS(PS-1.D0) > 1.D-7 ) write(0, *) "PS != 1: ", PS 
1745:  
1746:          ISTEP = ISTEP + 1 
1747:          ICNT  = ICNT + 1 
1748:  
1749: !     ADVANCE POSITIONS, ORIENTATIONS, AND THEIR TIME DERIVATIVES 
1750: !     MOVE also calls CHECK() and FORQUE 
1751:          
1752:          CALL MOVE(ISTEP, NEQ, PE, VRL, TKE, RKE) 
1753:  
1754: !     CALCULATE TEMPERATURES, PRESSURE AND TOTEL ENERGY AT TIME t 
1755:  
1756:          TMPTR   = 2.D0 * TKE  / DBLE(3 * NMOL - 3) 
1757:          TMPRT   = 2.D0 * RKE / DBLE(3 * NMOL)            
1758:          KE      = TKE + RKE 
1759:          !DH is used only in NVT 
1760:          DH      = ABS((KE + PE + 0.5D0*PS*PS/MQ + G*TMPFIX*LOG(S) - HNOT)/HNOT)  
1761:          TMP     = 2.D0 * KE / DBLE(G) 
1762:          PEPP    = PE / DBLE(NMOL) 
1763:          KEPP    = KE / DBLE(NMOL) 
1764:          EPP     = PEPP + KEPP 
1765:          PRS     = (2.D0 * TKE + VRL) / (3.D0 * VLM) 
1766:          SUME    = SUME + EPP 
1767:          SUMKE   = SUMKE + KEPP 
1768:          SUMPE   = SUMPE + PEPP 
1769:          SUMTMT  = SUMTMT + TMPTR 
1770:          SUMTMR  = SUMTMR + TMPRT 
1771:          SUMTMP  = SUMTMP + TMP 
1772:          SUMPRS  = SUMPRS + PRS 
1773:  
1774:          !do block averaging of potential energy, temperature, and pressure 
1775:          IF (ISTEP > NEQ) THEN 
1776:  
1777:             BSUMPE  = BSUMPE + PEPP 
1778:             BSUMTMP = BSUMTMP + TMP 
1779:             BSUMPRS = BSUMPRS + PRS 
1780:  
1781:             IF (MOD((ISTEP-NEQ),IDUMP4) == 0) THEN 
1782:  
1783:                AVPEB  = BSUMPE / DBLE(IDUMP4) 
1784:                AVTMPB = BSUMTMP / DBLE(IDUMP4) 
1785:                AVPRSB = BSUMPRS / DBLE(IDUMP4)   
1786:  
1787:                OPEN (UNIT=25, FILE='blockav1.dat', STATUS='UNKNOWN', POSITION ='APPEND') 
1788:                  WRITE (25, *) AVPEB, AVTMPB, AVPRSB 
1789:                CLOSE (UNIT=25, STATUS='KEEP') 
1790:  
1791:                BSUMPE  = 0.D0 
1792:                BSUMTMP = 0.D0 
1793:                BSUMPRS = 0.D0 
1794:  
1795:             ENDIF 
1796:  
1797:          ENDIF 
1798:   
1799:          !write out energies, temperatures, and their averages, etc. 
1800:          IF (MOD(ISTEP, IDUMP1) .EQ. 0) THEN 
1801:  
1802:             AVE    = SUME / DBLE(ICNT) 
1803:             AVKE   = SUMKE / DBLE(ICNT) 
1804:             AVPE   = SUMPE / DBLE(ICNT) 
1805:             AVTMPT = SUMTMT / DBLE(ICNT) 
1806:             AVTMPR = SUMTMR / DBLE(ICNT) 
1807:             AVTMP  = SUMTMP / DBLE(ICNT) 
1808:             AVPRS  = SUMPRS / DBLE(ICNT) 
1809:  
1810:             OPEN (UNIT=3, FILE='enrg1.dat', STATUS='UNKNOWN', POSITION ='APPEND') 
1811:                  WRITE (3, *) KEPP, PEPP, EPP  
1812:             CLOSE (UNIT=3, STATUS='KEEP') 
1813:             OPEN (UNIT=4, FILE='tmpprs1.dat', STATUS = 'UNKNOWN', POSITION ='APPEND') 
1814:                  WRITE (4, *) TMP, PRS 
1815:             CLOSE (UNIT=4, STATUS='KEEP')  
1816:  
1817:             OPEN (UNIT=33, FILE='avenrg1.dat', STATUS='UNKNOWN', POSITION ='APPEND') 
1818:                  WRITE (33, *) AVKE, AVPE, AVE 
1819:             CLOSE (UNIT=33, STATUS='KEEP') 
1820:             OPEN (UNIT=34, FILE='avtmpprs1.dat', STATUS = 'UNKNOWN', POSITION ='APPEND') 
1821:                  WRITE (34, *) AVTMP, AVPRS 
1822:             CLOSE (UNIT=34, STATUS='KEEP')  
1823:  
1824:             OPEN (UNIT=35, FILE='avtmptr1.dat', STATUS = 'UNKNOWN', POSITION ='APPEND') 
1825:                  WRITE (35, *) AVTMPT, AVTMPR 
1826:             CLOSE (UNIT=35, STATUS='KEEP') 
1827:  
1828:             OPEN (UNIT=36, FILE='error1.dat', STATUS = 'UNKNOWN', POSITION ='APPEND') 
1829:                  WRITE (36, *) DH !DH is NVT only 
1830:             CLOSE (UNIT=36, STATUS='KEEP') 
1831:  
1832:          ENDIF 
1833:  
1834:          ! calculate energy averaged over time and molecules. 
1835:          ! calculate variance over molecules of time averaged energy 
1836:          IF (ISTEP > NEQ) THEN 
1837:  
1838:             CALL EFM (EFMT, TIME) 
1839:  
1840:          ENDIF 
1841:  
1842:          !scale velocities. NVE only 
1843:          IF (ISTEP <= NBATH .AND. MOD(ISTEP, IDUMP3) == 0 .AND. ISTEP>0) THEN 
1844:  
1845:             CALL NVESCLVEL(AVTMPT, AVTMPR) 
1846:             ICNT   = 0 
1847:             SUME   = 0.D0 
1848:             SUMKE  = 0.D0 
1849:             SUMPE  = 0.D0 
1850:             SUMTMP = 0.D0 
1851:             SUMTMT = 0.D0 
1852:             SUMTMR = 0.D0 
1853:             SUMPRS = 0.D0 
1854:  
1855:          ENDIF 
1856:  
1857:          !write out positions, orientations, and angular change since t=NEQ 
1858:          IF (ISTEP > NEQ .AND. MOD(ISTEP, IDUMP2) == 0) THEN    
1859:  
1860:            OPEN (UNIT=7,  FILE='pos1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
1861:            OPEN (UNIT=8,  FILE='ortn1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
1862:            OPEN (UNIT=28, FILE='rot1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
1863:            OPEN (UNIT=29, FILE='efm1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
1864:  
1865:            DO I = 1, NMOL 
1866:  
1867:              WRITE (7, *)  R(I,1), R(I,2), R(I,3) 
1868:              WRITE (8, 900)  QTRN(I,1), QTRN(I,2), QTRN(I,3), QTRN(I,4)  
1869:              WRITE (28, *) PHI(I,1), PHI(I,2), PHI(I,3)   
1870:              900 FORMAT (1x,4ES25.16) 
1871:  
1872:            ENDDO 
1873:  
1874:            WRITE(29, *) TIME, EFMT 
1875:  
1876:            CLOSE (UNIT=7,  STATUS='KEEP')   
1877:            CLOSE (UNIT=8,  STATUS='KEEP')   
1878:            CLOSE (UNIT=28, STATUS='KEEP') 
1879:            CLOSE (UNIT=29, STATUS='KEEP') 
1880:  
1881:            IF (MOD(ISTEP, IDUMP4) .EQ. 0) THEN    
1882:  
1883:              OPEN (UNIT=9, FILE='lvel1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
1884:              OPEN (UNIT=10, FILE='avel1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
1885:  
1886:              DO I = 1, NMOL 
1887:  
1888:                WRITE (9, *)  V(I,1)/S, V(I,2)/S, V(I,3)/S 
1889:                WRITE (10, *) W(I,1), W(I,2), W(I,3) 
1890:  
1891:              ENDDO 
1892:  
1893:              CLOSE (UNIT = 9, STATUS = 'KEEP')   
1894:              CLOSE (UNIT = 10, STATUS ='KEEP')   
1895:  
1896:            ENDIF 
1897:  
1898:          ENDIF 
1899:                 
1900:          !reinitialize energy averages 
1901:          IF (ISTEP == NBATH .OR. ISTEP == NEQ) THEN 
1902:                
1903:             !CALL NVESCLVEL(AVTMPT, AVTMPR) 
1904:             ICNT   = 0 
1905:             SUME   = 0.D0 
1906:             SUMKE  = 0.D0 
1907:             SUMPE  = 0.D0 
1908:             SUMTMT = 0.D0 
1909:             SUMTMR = 0.D0 
1910:             SUMTMP = 0.D0 
1911:             SUMPRS = 0.D0         
1912:             TIME   = 0.D0 
1913:  
1914:             BSUMPE  = 0.D0 
1915:             BSUMTMP = 0.D0 
1916:             BSUMPRS = 0.D0 
1917:                           
1918:          ENDIF   
1919:  
1920:       ENDDO 
1921:  
1922:       !program ended, print out the final data 
1923:  
1924:       OPEN (UNIT = 21, FILE = 'finalpos1.dat', STATUS = 'UNKNOWN') 
1925:       OPEN (UNIT = 22, FILE = 'finalortn1.dat', STATUS = 'UNKNOWN') 
1926:       OPEN (UNIT = 23, FILE = 'finallinvel1.dat', STATUS = 'UNKNOWN') 
1927:       OPEN (UNIT = 24, FILE = 'finalangvel1.dat', STATUS = 'UNKNOWN') 
1928:        
1929:       DO I = 1, NMOL 
1930:  
1931:          WRITE (21, *) R(I,1), R(I,2), R(I,3) 
1932:          WRITE (22, 900) QTRN(I,1), QTRN(I,2), QTRN(I,3), QTRN(I,4) 
1933:          WRITE (23, *) V(I,1), V(I,2), V(I,3) 
1934:          WRITE (24, *) W(I,1), W(I,2), W(I,3) 
1935:  
1936:       ENDDO 
1937:  
1938:       CLOSE(UNIT=21) 
1939:       CLOSE(UNIT=22) 
1940:       CLOSE(UNIT=23) 
1941:       CLOSE(UNIT=24) 
1942:  
1943: !      DEALLOCATE (R, QTRN, V, W, F, T, P, T4, PHI, FSITE, MST, KEP, PEP, SUMEPT, STAT = STATUS) 
1944:  
1945:       call RUN_INFO() 
1946:  
1947: !      END PROGRAM LWOTP_MD_NVE 
1948:       END SUBROUTINE RIGIDMD_NVE 
1949: !     ------------------------------------------------------------------------------------ 
1950:   
1951:       SUBROUTINE NVESTEP1 (S, PS, FCTMQT ) 
1952:  
1953:       IMPLICIT NONE 
1954:  
1955:       DOUBLE PRECISION, INTENT(IN) :: S, PS 
1956:       DOUBLE PRECISION, INTENT(IN) :: FCTMQT 
1957:  
1958:       END SUBROUTINE NVESTEP1 
1959:  
1960: !     ------------------------------------------------------------------------------------ 
1961:   
1962:       SUBROUTINE NVESTEP2 (PE) 
1963:  
1964:       !STEP 2 in the MOVE algorithm 
1965:       ! STEP 2: update V, P, PS 
1966:       ! PS = PS - PE*HALFDT   !! PE=potential energy 
1967:       ! loop over molecules 
1968:       !   update linear momentum (or velocity) 
1969:       !   V(I,:) = V(I,:) + S * F(I,:) * HALFDT / MASS    
1970:       !   update rotational momentum  
1971:       !   P(I,:) = P(I,:) +  2*S * DELT/2 * MATMUL(MX,T4(I,:)) 
1972:       ! end loop over molecules 
1973:        
1974:       USE MDCOMMONS, only : NMOL, HALFDT, MASS, QTRN, DELT, T4, P, V, F 
1975:  
1976:       IMPLICIT NONE 
1977:  
1978:       INTEGER :: I 
1979:       DOUBLE PRECISION :: MX(4,4), Q(4) 
1980:       DOUBLE PRECISION, INTENT(IN) :: PE 
1981:  
1982:       !PS     = PS - PE * HALFDT 
1983:  
1984:       DO I = 1, NMOL 
1985:  
1986:          V(I,:) = V(I,:) + F(I,:) * HALFDT / MASS           
1987:  
1988:          Q      = QTRN(I,:) 
1989:  
1990:          CALL CONSTRUCT_MX( MX, Q ) 
1991:  
1992:          P(I,:)     = P(I,:) +  DELT * MATMUL(MX,T4(I,:)) 
1993:  
1994:       ENDDO 
1995:  
1996:       END SUBROUTINE NVESTEP2 
1997:  
1998: !     ------------------------------------------------------------------------------------ 
1999:   
2000:       SUBROUTINE NVESTEP3 () 
2001:  
2002:       !STEP 3 in the MOVE algorithm 
2003:       ! STEP 3: update QTRN, P, PS with contributions from principle axis 3 
2004:       ! ZETASUM = 0 
2005:       ! loop over molecules 
2006:       !   ZETA(I,3) = dot_product( P(I,:) , MX(:,4) ) / ( 4 * S * JMI(3) ) 
2007:       !   ZETASUM = ZETASUM + ZETA(I,3)**2 
2008:       !   Q(I) = cos( ZETA(I,3)*HALFDT )*Q(I) + sin( ZETA(I,3)*HALFDT )*MX(:,4) 
2009:       !   P(I) = cos( ZETA(I,3)*HALFDT )*P(I) + sin( ZETA(I,3)*HALFDT )*MP(:,4) 
2010:       !   !   in the above, MP(:,:) is the analogous matrix to MX for vector P(I,:), 
2011:       !   !   i.e. MP(:,4) = (/-P(I,4),P(I,3),-P(I,2),P(I,1)/) 
2012:       ! end loop over molecules 
2013:       ! PS = PS + ZETASUM*JMI(3)*DELT 
2014:       ! contribution from principle axis 2 
2015:        
2016:       USE MDCOMMONS, only : QTRN, P, HALFDT, JMI, NMOL  
2017:  
2018:       IMPLICIT NONE 
2019:  
2020:       INTEGER :: I 
2021:       DOUBLE PRECISION :: SMPHSQ, Q(4), DQ(4), PHIT 
2022:       DOUBLE PRECISION :: PI(4), CP, SP 
2023:  
2024:       SMPHSQ = 0.D0 
2025:       DO I = 1, NMOL 
2026:  
2027:          Q         = QTRN(I,:) 
2028:          PI        = P(I,:) 
2029:          !DQ        = MX(:,4) 
2030:          DQ        = (/-Q(4),Q(3),-Q(2),Q(1)/) 
2031:          PHIT      = 0.25D0*DOT_PRODUCT(PI,DQ)*HALFDT/(JMI(3)) 
2032:          !SMPHSQ    = SMPHSQ + PHIT * PHIT 
2033:          CP        = COS(PHIT) 
2034:          SP        = SIN(PHIT) 
2035:          P(I,:)    = CP*PI + SP*(/-PI(4),PI(3),-PI(2),PI(1)/) 
2036:          QTRN(I,:) = CP*Q + SP*DQ 
2037:  
2038:       ENDDO 
2039:  
2040:       !PS     = PS + JMI(3)*SMPHSQ*DELT 
2041:  
2042:       END SUBROUTINE NVESTEP3 
2043:  
2044: !     ------------------------------------------------------------------------------------ 
2045:   
2046:       SUBROUTINE NVESTEP4 () 
2047:  
2048:       !STEP 4 in the MOVE algorithm 
2049:       ! STEP 4: update QTRN, P, PS with contributions from principle axis 2 
2050:       ! ZETASUM = 0 
2051:       ! loop over molecules 
2052:       !   ZETA(I,2) = dot_product( P(I,:) , MX(:,4) ) / ( 4 * S * JMI(2) ) 
2053:       !   ZETASUM = ZETASUM + ZETA(I,2)**2 
2054:       !   Q(I) = cos( ZETA(I,2)*HALFDT )*Q(I) + sin( ZETA(I,2)*HALFDT )*MX(:,2) 
2055:       !   P(I) = cos( ZETA(I,2)*HALFDT )*P(I) + sin( ZETA(I,2)*HALFDT )*MP(:,2) 
2056:       !   !   in the above, MP(:,:) is the analogous matrix to MX for vector P(I,:), 
2057:       !   !   i.e. MP(:,2) = (/-P(I,3),-P(I,4),P(I,1),P(I,2)/) 
2058:       ! end loop over molecules 
2059:       ! PS = PS + ZETASUM*JMI(2)*DELT 
2060:        
2061:       USE MDCOMMONS, only : QTRN, P, HALFDT, JMI, NMOL  
2062:  
2063:       IMPLICIT NONE 
2064:  
2065:       INTEGER :: I 
2066:       DOUBLE PRECISION :: SMPHSQ, Q(4), DQ(4), PHIT 
2067:       DOUBLE PRECISION :: PI(4), CP, SP 
2068:  
2069:       SMPHSQ = 0.D0 
2070:       DO I = 1, NMOL 
2071:  
2072:          Q         = QTRN(I,:) 
2073:          PI        = P(I,:) 
2074:          DQ        = (/-Q(3),-Q(4),Q(1),Q(2)/) 
2075:          PHIT      = 0.25D0*DOT_PRODUCT(PI,DQ)*HALFDT/(JMI(2)) 
2076:          SMPHSQ    = SMPHSQ + PHIT*PHIT 
2077:          CP        = COS(PHIT) 
2078:          SP        = SIN(PHIT) 
2079:          P(I,:)    = CP*PI + SP*(/-PI(3),-PI(4),PI(1),PI(2)/) 
2080:          QTRN(I,:) = CP*Q + SP*DQ 
2081:  
2082:       ENDDO 
2083:  
2084:       !PS     = PS + JMI(2)*SMPHSQ*DELT 
2085:  
2086:       END SUBROUTINE NVESTEP4 
2087:  
2088: !     ------------------------------------------------------------------------------------ 
2089:   
2090:       SUBROUTINE NVESTEP5 () 
2091:  
2092:       !STEP 5 in the MOVE algorithm 
2093:       ! STEP 5: update QTRN, P, PS with contributions from principle axis 1 
2094:       ! STEP 5 is analogous to STEP 3 and STEP 4 but for ZETA(:,1), i.e. the 
2095:       ! first principle axis, 
2096:       ! EXCEPT, the update of PS is different and involves the translational 
2097:       ! kinetic energy, the temperature (TMPFIX), the number of degrees of 
2098:       ! freedom (g), and HNOT, 
2099:        
2100:       USE MDCOMMONS, only : QTRN, P, HALFDT, JMI, NMOL, DELT, R, V 
2101:  
2102:       IMPLICIT NONE 
2103:  
2104:       INTEGER :: I 
2105:       DOUBLE PRECISION :: SMPHSQ, Q(4), DQ(4), PHIT 
2106:       DOUBLE PRECISION :: PI(4), CP, SP 
2107:       !DOUBLE PRECISION :: TKE !translational kinetic energy 
2108:  
2109:       SMPHSQ = 0.D0 
2110:       !TKE = 0 
2111:       DO I = 1, NMOL 
2112:  
2113:          R(I,:)    = R(I,:) + V(I,:) * DELT 
2114:          Q         = QTRN(I,:) 
2115:          PI        = P(I,:) 
2116:          DQ        = (/-Q(2),Q(1),Q(4),-Q(3)/) 
2117:          PHIT      = 0.25D0*DOT_PRODUCT(PI,DQ)*DELT/(JMI(1)) 
2118:          !SMPHSQ    = SMPHSQ + PHIT * PHIT 
2119:          CP        = COS(PHIT) 
2120:          SP        = SIN(PHIT) 
2121:          P(I,:)    = CP*PI + SP*(/-PI(2),PI(1),PI(4),-PI(3)/) 
2122:          QTRN(I,:) = CP*Q  + SP*DQ 
2123:          !TKE       = TKE + V(I,1)*V(I,1) + V(I,2)*V(I,2) + V(I,3)*V(I,3) 
2124:  
2125:       ENDDO 
2126:  
2127:       !PS = PS + (0.5D0*MASS*TKE/(S*S) + 2.D0*JMI(1)*SMPHSQ - DBLE(G)*TMPFIX*LOG(S) & 
2128:          !+ HNOT - DBLE(G)*TMPFIX) * DELT 
2129:  
2130:       END SUBROUTINE NVESTEP5 
2131:  
2132: !     ---------------------------------------------------------------------------------------------- 
2133:  
2134:       SUBROUTINE NVESCLVEL(TMPT, TMPR) 
2135:  
2136: ! scale linear velocities??? 
2137:  
2138:       USE MDCOMMONS, only : TMPFIX, V, W, QTRN, NMOL, JMI, P 
2139:  
2140:       IMPLICIT NONE 
2141:  
2142:       INTEGER          :: I 
2143:       DOUBLE PRECISION :: SCLTR, SCLRT, TMPT, TMPR 
2144:       DOUBLE PRECISION :: Q(4), RMXT(3,3), MX(4,4), W4(4) 
2145:  
2146:       SCLTR   = DSQRT(TMPFIX / TMPT) 
2147:       SCLRT   = DSQRT(TMPFIX / TMPR) 
2148:  
2149:       DO I = 1, NMOL 
2150:  
2151:          V(I,:) = V(I,:) * SCLTR 
2152:  
2153:          W(I,:) = W(I,:) * SCLRT 
2154:  
2155:          Q    = QTRN(I,:) 
2156:          CALL CONSTRUCT_RMXT(RMXT,Q) 
2157:  
2158: !     CONSTRUCT THE MATRIX MX FROM THE QUATERNION Q 
2159:  
2160:          CALL CONSTRUCT_MX(MX,Q) 
2161:  
2162:          W4(1)  = 0.D0 
2163:          W4(2)  = 2.D0*JMI(1)*W(I,1) 
2164:          W4(3)  = 2.D0*JMI(2)*W(I,2) 
2165:          W4(4)  = 2.D0*JMI(3)*W(I,3) 
2166:  
2167:          P(I,:) = MATMUL(MX,W4) 
2168:  
2169:       ENDDO 
2170:  
2171:       END SUBROUTINE NVESCLVEL 
2172: ! 
2173: ! svn info  
2174: ! $Date: 2011-09-21 11:24:09 +0100 (Wed, 21 Sep 2011) $ 
2175: ! $Rev: 25 $ 
2176: ! 
2177:  
2178: ! the lewis wahnstrom model has  
2179:       ! sigma = 0.483 nm 
2180:       ! epsilon = 5.276 kJ/mol 
2181:   
2182: !     ---------------------------------------------------------------------------------------------- 
2183:  
2184: !      PROGRAM LWOTP_MD_NVE 
2185:       SUBROUTINE RIGIDMD_NVE_CHANGE_DENSITY 
2186:       USE MDCOMMONS 
2187:  
2188:       IMPLICIT NONE 
2189:  
2190: ! IDUMP1: dump files.  always:  
2191: !                     enrg1.dat: KEPP, PEPP, EPP: kinetic, potential, energy per molecule 
2192: !                     tmpprs1.dat: temp, pressure? 
2193: !                     avenrg1.dat: average KE, PE, E 
2194: !                     avtmpprs1.dat: average temp, pressure? 
2195: !                     avtmptr1.dat: average tranlational, rotational temperature 
2196: ! IDUMP2: dump files.  ISTEP > NEQ 
2197: !                     pos1.dat:  xyz positions of CoM of molecules 
2198: !                     ortn1.dat: quaternions giving the orientation of molecules 
2199: !                     rot1.dat:  phi1, phi2, phi3 ???? 
2200: !                     efm1.dat:  energy fluctuation metric: time from NEQ, and variance over molecules of time averaged energy  
2201: ! IDUMP3: CALL SCLVEL.  ISTEP <= NBATH: scale linear velocities ??? 
2202: ! IDUMP4: dump files.  ISTEP > NEQ:  
2203: !                     blockav1.dat: average PE, temp, pressure 
2204: !                     lvel1.dat:    linear velocities 
2205: !                     avel1.dat:    angular velocities 
2206: ! TKE: translational kinetic energy 
2207: ! RKE: rotational kinetic energy 
2208: ! MX:  is a representation of the quaternion as a matrix in such a way 
2209: !      that quaternion addition and multiplication correspond to matrix 
2210: !      addition and matrix multiplication 
2211: !      In this way the conjugate of the quaternion corresponds to the 
2212: !      transpose of the matrix.  The fourth power of the norm of a 
2213: !      quaternion is the determinant of the corresponding matrix. 
2214: !      Complex numbers are block diagonal matrices with two 2x2 blocks. 
2215: ! EFMT: Energy Fluctuation Metric: the variance over molecules of the 
2216: !       time averaged energy 
2217:       INTEGER          :: NSTEP, NBATH, NEQ, ISTEP, ICNT, IDUMP1, IDUMP2, IDUMP3, IDUMP4, I!, K!, STATUS 
2218:       DOUBLE PRECISION :: RCUT6, RCUT12, TMP, TMPTR, TMPRT, RHO, VLM 
2219:       DOUBLE PRECISION :: Q(4), RMXT(3,3), MX(4,4), W4(4) 
2220:       DOUBLE PRECISION :: PRS, SUMPRS, AVPRS, PE, KE, TKE, RKE, PEPP, KEPP, EPP, VRL 
2221:       DOUBLE PRECISION :: SUME, SUMKE, SUMPE, SUMTMT, SUMTMR, SUMTMP 
2222:       DOUBLE PRECISION :: AVE, AVKE, AVPE, AVTMPT, AVTMPR, AVTMP, EFMT, TIME 
2223:       DOUBLE PRECISION :: BSUMPE, BSUMTMP, BSUMPRS, AVPEB, AVTMPB, AVPRSB 
2224:       DOUBLE PRECISION :: MQ=0.1D0, DH=1.D0 !only NVT 
2225:       INTEGER :: TMPINT !initial temperature 
2226:       logical :: bool1, bool2 
2227:       !CALL RANDOM_SEED() 
2228:       DOUBLE PRECISION :: RHOFINAL, DELRHO, RHOOLD, LCHANGE, DBLE 
2229:  
2230:       OPEN (UNIT = 1, FILE = 'parameter1.inp', STATUS = 'UNKNOWN') 
2231:  
2232: !     READ INPUT PARAMETERS 
2233:  
2234:       READ (1, *) 
2235:       READ (1, *) NMOL, NSITE 
2236:       READ (1, *)  
2237:       READ (1, *) RCUT, RLIST 
2238:       READ (1, *) 
2239:       READ (1, *) TMPFIX 
2240:       READ (1, *) 
2241:       READ (1, *) RHO, RHOFINAL, DELRHO 
2242:       READ (1, *) 
2243:       READ (1, *) DELT 
2244:       READ (1, *) 
2245:       READ (1, *) NSTEP, NBATH, NEQ 
2246:       READ (1, *)  
2247:       READ (1, *) IDUMP1, IDUMP2, IDUMP3, IDUMP4 
2248:   
2249:       CLOSE (1) 
2250:       TMPINT = TMPFIX !initial temperature 
2251:  
2252: !     CALCULATE FROM INPUT PARAMETERS 
2253:  
2254:       G      = 6 * NMOL - 3 
2255:       NTST   = NMOL * NSITE             ! TOTAL NUMBER OF SITES 
2256:       RCUTSQ = RCUT * RCUT 
2257:       RLSTSQ = RLIST * RLIST 
2258:       EPS4   = 4.D0  ! 4*epsilon 
2259:       RCUT6  = (1.D0 / RCUT) ** 6 
2260:       RCUT12 = RCUT6 * RCUT6 
2261:       CNSTA  = (6.D0 * RCUT12 - 3.D0 * RCUT6) / RCUTSQ 
2262:       CNSTB  = 4.D0 * RCUT6 - 7.D0 * RCUT12 
2263:       HALFDT = 0.5D0 * DELT 
2264:       FCTMQT = 0.5D0*HALFDT/MQ !used only in NVT 
2265:   
2266: !     ALLOCATE THE ARRAYS 
2267:  
2268: !      ALLOCATE (R(NMOL,3), QTRN(NMOL,4), V(NMOL,3), W(NMOL,3), F(NMOL,3), T(NMOL,3), P(NMOL,4), & 
2269: !      T4(NMOL,4), PHI(NMOL,3), FSITE(NTST,3), MST(NSITE,3), KEP(NMOL), PEP(NMOL), SUMEPT(NMOL), & 
2270: !      STAT = STATUS) 
2271:  
2272: !      IF (STATUS /= 0) STOP 'ALLOCATION PROBLEM' 
2273:   
2274: !     BOX LENGTH 
2275:  
2276:       VLM    = DBLE(NMOL) / RHO 
2277:       BOXL   = VLM ** (1.D0/3.D0) 
2278:  
2279:       CALL DEFMOL() 
2280:  
2281: !      WRITE(*,*) RCUT, RCUT2, RLIST 
2282:  
2283:  
2284:       INQUIRE(FILE="initpos1.dat", EXIST=bool1 )  
2285:       INQUIRE(FILE="initortn1.dat", EXIST=bool2 )  
2286:       if ( bool1 .and. bool2 ) then 
2287:         !start from an existing configuration 
2288:         !CALL EQCON() 
2289:         write(*,*) " reading init files" 
2290:         CALL LOAD_POS_ORTN() 
2291:       else 
2292:         !     start from scratch 
2293:         !     START WITH A SIMPLE CUBIC LATTICE CONFIGURATION 
2294:         write(*,*) " no init files: starting from scratch" 
2295:         CALL RIGIDMD_SC() 
2296:         !     SET INITIAL ORIENTATIONS 
2297:         CALL INTORN() 
2298:       endif 
2299:  
2300: !     SET INITIAL LINEAR and angular VELICITIES 
2301:       CALL INTVEL(TMPINT) 
2302:  
2303:  
2304: !     START SIMULATION 
2305:  
2306:       ISTEP  = 0   
2307:       ICNT   = 0 
2308:       SUME   = 0.D0 
2309:       SUMKE  = 0.D0 
2310:       SUMPE  = 0.D0 
2311:       SUMTMT = 0.D0 
2312:       SUMTMR = 0.D0 
2313:       SUMTMP = 0.D0 
2314:       SUMPRS = 0.D0 
2315:        
2316:       PHI(:,:) = 0.D0 
2317:       KE       = 0.D0 
2318:  
2319: !     CALCULATE FORCES AND TORQUES AT TIME t=0 
2320:  
2321:       UPDATE = .TRUE. 
2322:  
2323:       CALL FORQUE(PE, VRL) 
2324:  
2325:       DO I = 1, NMOL 
2326:  
2327:          !calculate the transposed rotation matrix from quaternion 
2328:          Q    = QTRN(I,:) 
2329:          CALL CONSTRUCT_RMXT( RMXT, Q ) 
2330:  
2331: !        TORQUE AT TIME t=0 IN THE BODY FRAME USING THE TRANSPOSE OF THE ROTATION MATRIX 
2332:  
2333:          T4(I,1)   = 0.D0 
2334:          T4(I,2:4) = MATMUL(RMXT,T(I,:)) 
2335:  
2336: !     CONSTRUCT THE MATRIX MX FROM THE QUATERNION Q 
2337:  
2338:          CALL CONSTRUCT_MX( MX, Q ) 
2339:  
2340:          W4(1)  = 0.D0 
2341:          W4(2)  = 2.D0*W(I,1) 
2342:          W4(3)  = 2.D0*W(I,2) 
2343:          W4(4)  = 2.D0*W(I,3) 
2344:  
2345:          ! intitialize the momentum conjugate to angular velocity 
2346:          P(I,:) = MATMUL(MX,W4) 
2347:   
2348:          KE = KE + 0.5D0*MASS*(V(I,1)*V(I,1) + V(I,2)*V(I,2) + V(I,3)*V(I,3)) & 
2349:                  + 0.5D0*(W(I,1)*W(I,1)/JMI(1) + W(I,2)*W(I,2)/JMI(2) + W(I,3)*W(I,3)/JMI(3)) 
2350:  
2351:       ENDDO 
2352:  
2353:       HNOT = KE + PE  
2354:  
2355:       DO WHILE (ISTEP < NSTEP)             
2356:  
2357:          !check to make sure the extra variables relevant only for NVT runs are 
2358:          !constant and = 1 
2359:          !UNIT 0 should be standard error 
2360:          IF ( ABS(S-1.D0) > 1.D-7 ) write(0, *) "S != 1: ", S 
2361:          IF ( ABS(PS-1.D0) > 1.D-7 ) write(0, *) "PS != 1: ", PS 
2362:  
2363:          ISTEP = ISTEP + 1 
2364:          ICNT  = ICNT + 1 
2365:  
2366: !     ADVANCE POSITIONS, ORIENTATIONS, AND THEIR TIME DERIVATIVES 
2367: !     MOVE also calls CHECK() and FORQUE 
2368:          
2369:          CALL MOVE(ISTEP, NEQ, PE, VRL, TKE, RKE) 
2370:  
2371: !     CALCULATE TEMPERATURES, PRESSURE AND TOTEL ENERGY AT TIME t 
2372:  
2373:          TMPTR   = 2.D0 * TKE  / DBLE(3 * NMOL - 3) 
2374:          TMPRT   = 2.D0 * RKE / DBLE(3 * NMOL)            
2375:          KE      = TKE + RKE 
2376:          !DH is used only in NVT 
2377:          DH      = ABS((KE + PE + 0.5D0*PS*PS/MQ + G*TMPFIX*LOG(S) - HNOT)/HNOT)  
2378:          TMP     = 2.D0 * KE / DBLE(G) 
2379:          PEPP    = PE / DBLE(NMOL) 
2380:          KEPP    = KE / DBLE(NMOL) 
2381:          EPP     = PEPP + KEPP 
2382:          PRS     = (2.D0 * TKE + VRL) / (3.D0 * VLM) 
2383:          SUME    = SUME + EPP 
2384:          SUMKE   = SUMKE + KEPP 
2385:          SUMPE   = SUMPE + PEPP 
2386:          SUMTMT  = SUMTMT + TMPTR 
2387:          SUMTMR  = SUMTMR + TMPRT 
2388:          SUMTMP  = SUMTMP + TMP 
2389:          SUMPRS  = SUMPRS + PRS 
2390:  
2391:          !do block averaging of potential energy, temperature, and pressure 
2392:          IF (ISTEP > NEQ) THEN 
2393:  
2394:             BSUMPE  = BSUMPE + PEPP 
2395:             BSUMTMP = BSUMTMP + TMP 
2396:             BSUMPRS = BSUMPRS + PRS 
2397:  
2398:             IF (MOD((ISTEP-NEQ),IDUMP4) == 0) THEN 
2399:  
2400:                AVPEB  = BSUMPE / DBLE(IDUMP4) 
2401:                AVTMPB = BSUMTMP / DBLE(IDUMP4) 
2402:                AVPRSB = BSUMPRS / DBLE(IDUMP4)   
2403:  
2404:                OPEN (UNIT=25, FILE='blockav1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
2405:                  WRITE (25, *) AVPEB, AVTMPB, AVPRSB 
2406:                CLOSE (UNIT=25, STATUS='KEEP') 
2407:  
2408:                BSUMPE  = 0.D0 
2409:                BSUMTMP = 0.D0 
2410:                BSUMPRS = 0.D0 
2411:  
2412:             ENDIF 
2413:  
2414:          ENDIF 
2415:   
2416:          !write out energies, temperatures, and their averages, etc. 
2417:          IF (MOD(ISTEP, IDUMP1) .EQ. 0) THEN 
2418:  
2419:             AVE    = SUME / DBLE(ICNT) 
2420:             AVKE   = SUMKE / DBLE(ICNT) 
2421:             AVPE   = SUMPE / DBLE(ICNT) 
2422:             AVTMPT = SUMTMT / DBLE(ICNT) 
2423:             AVTMPR = SUMTMR / DBLE(ICNT) 
2424:             AVTMP  = SUMTMP / DBLE(ICNT) 
2425:             AVPRS  = SUMPRS / DBLE(ICNT) 
2426:  
2427:             OPEN (UNIT=3, FILE='enrg1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
2428:                  WRITE (3, *) KEPP, PEPP, EPP  
2429:             CLOSE (UNIT=3, STATUS='KEEP') 
2430:             OPEN (UNIT=4, FILE='tmpprs1.dat', STATUS = 'UNKNOWN', POSITION='APPEND') 
2431:                  WRITE (4, *) TMP, PRS 
2432:             CLOSE (UNIT=4, STATUS='KEEP')  
2433:  
2434:             OPEN (UNIT=33, FILE='avenrg1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
2435:                  WRITE (33, *) AVKE, AVPE, AVE 
2436:             CLOSE (UNIT=33, STATUS='KEEP') 
2437:             OPEN (UNIT=34, FILE='avtmpprs1.dat', STATUS = 'UNKNOWN', POSITION='APPEND') 
2438:                  WRITE (34, *) AVTMP, AVPRS 
2439:             CLOSE (UNIT=34, STATUS='KEEP')  
2440:  
2441:             OPEN (UNIT=35, FILE='avtmptr1.dat', STATUS = 'UNKNOWN', POSITION='APPEND') 
2442:                  WRITE (35, *) AVTMPT, AVTMPR 
2443:             CLOSE (UNIT=35, STATUS='KEEP') 
2444:  
2445:             OPEN (UNIT=36, FILE='error1.dat', STATUS = 'UNKNOWN', POSITION='APPEND') 
2446:                  WRITE (36, *) DH !DH is NVT only 
2447:             CLOSE (UNIT=36, STATUS='KEEP') 
2448:  
2449:          ENDIF 
2450:  
2451:          ! calculate energy averaged over time and molecules. 
2452:          ! calculate variance over molecules of time averaged energy 
2453:          IF (ISTEP > NEQ) THEN 
2454:  
2455:             CALL EFM (EFMT, TIME) 
2456:  
2457:          ENDIF 
2458:  
2459:          !scale velocities. NVE only 
2460:          IF ( MOD(ISTEP, IDUMP3) == 0) THEN 
2461:  
2462:              RHOOLD=RHO 
2463:            IF ( RHOFINAL .gt. RHO+DELRHO ) THEN 
2464:              RHO=RHO+DELRHO 
2465:            ELSE IF ( RHOFINAL .lt. RHO-DELRHO ) THEN 
2466:              RHO=RHO-DELRHO 
2467:            ELSE 
2468:              RHO=RHOFINAL 
2469:            ENDIF 
2470:            if ( RHOOLD .ne. RHO ) then 
2471:              VLM    = DBLE(NMOL) / RHO 
2472:              BOXL   = VLM ** (1.D0/3.D0) 
2473:  
2474:              LCHANGE=(RHOOLD/RHO)**(1.D0/3.D0) 
2475:              write (*,*) "RHO ", rho, " LCHANGE ", LCHANGE 
2476:              R = R*LCHANGE; 
2477:            endif 
2478:  
2479:  
2480:             CALL NVESCLVEL(AVTMPT, AVTMPR) 
2481:             ICNT   = 0 
2482:             SUME   = 0.D0 
2483:             SUMKE  = 0.D0 
2484:             SUMPE  = 0.D0 
2485:             SUMTMP = 0.D0 
2486:             SUMTMT = 0.D0 
2487:             SUMTMR = 0.D0 
2488:             SUMPRS = 0.D0 
2489:  
2490:          ENDIF 
2491:  
2492:          !write out positions, orientations, and angular change since t=NEQ 
2493:          IF (ISTEP > NEQ .AND. MOD(ISTEP, IDUMP2) == 0) THEN    
2494:  
2495:            OPEN (UNIT=7,  FILE='pos1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
2496:            OPEN (UNIT=8,  FILE='ortn1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
2497:            !OPEN (UNIT=28, FILE='rot1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
2498:            OPEN (UNIT=29, FILE='efm1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
2499:  
2500:            DO I = 1, NMOL 
2501:  
2502:              WRITE (7, *)  R(I,1), R(I,2), R(I,3) 
2503:              WRITE (8, 900)  QTRN(I,1), QTRN(I,2), QTRN(I,3), QTRN(I,4)  
2504:              !WRITE (28, *) PHI(I,1), PHI(I,2), PHI(I,3)   
2505:              900 FORMAT (1x,ES20.12,1x,ES20.12,1x,ES20.12) 
2506:  
2507:            ENDDO 
2508:  
2509:            WRITE(29, *) TIME, EFMT 
2510:  
2511:            CLOSE (UNIT=7,  STATUS='KEEP')   
2512:            CLOSE (UNIT=8,  STATUS='KEEP')   
2513:            !CLOSE (UNIT=28, STATUS='KEEP') 
2514:            CLOSE (UNIT=29, STATUS='KEEP') 
2515:  
2516:            IF (MOD(ISTEP, IDUMP4) .EQ. 0) THEN    
2517:  
2518:              OPEN (UNIT=9, FILE='lvel1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
2519:              OPEN (UNIT=10, FILE='avel1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
2520:  
2521:              DO I = 1, NMOL 
2522:  
2523:                WRITE (9, *)  V(I,1)/S, V(1,2)/S, V(1,3)/S 
2524:                WRITE (10, *) W(I,1), W(I,2), W(I,3) 
2525:  
2526:              ENDDO 
2527:  
2528:              CLOSE (UNIT = 9, STATUS = 'KEEP')   
2529:              CLOSE (UNIT = 10, STATUS ='KEEP')   
2530:  
2531:            ENDIF 
2532:  
2533:          ENDIF 
2534:                 
2535:          !reinitialize energy averages 
2536:          IF (ISTEP == NBATH .OR. ISTEP == NEQ) THEN 
2537:                
2538:             ICNT   = 0 
2539:             SUME   = 0.D0 
2540:             SUMKE  = 0.D0 
2541:             SUMPE  = 0.D0 
2542:             SUMTMT = 0.D0 
2543:             SUMTMR = 0.D0 
2544:             SUMTMP = 0.D0 
2545:             SUMPRS = 0.D0         
2546:             TIME   = 0.D0 
2547:  
2548:             BSUMPE  = 0.D0 
2549:             BSUMTMP = 0.D0 
2550:             BSUMPRS = 0.D0 
2551:                           
2552:          ENDIF   
2553:  
2554:       ENDDO 
2555:  
2556:       !program ended, print out the final data 
2557:  
2558:       OPEN (UNIT = 21, FILE = 'finalpos1.dat', STATUS = 'UNKNOWN') 
2559:       OPEN (UNIT = 22, FILE = 'finalortn1.dat', STATUS = 'UNKNOWN') 
2560:       OPEN (UNIT = 23, FILE = 'finallinvel1.dat', STATUS = 'UNKNOWN') 
2561:       OPEN (UNIT = 24, FILE = 'finalangvel1.dat', STATUS = 'UNKNOWN') 
2562:        
2563:       DO I = 1, NMOL 
2564:  
2565:          WRITE (21, *) R(I,1), R(I,2), R(I,3) 
2566:          WRITE (22, 900) QTRN(I,1), QTRN(I,2), QTRN(I,3), QTRN(I,4) 
2567:          WRITE (23, *) V(I,1), V(I,2), V(I,3) 
2568:          WRITE (24, *) W(I,1), W(I,2), W(I,3) 
2569:  
2570:       ENDDO 
2571:  
2572:       CLOSE(UNIT=21) 
2573:       CLOSE(UNIT=22) 
2574:       CLOSE(UNIT=23) 
2575:       CLOSE(UNIT=24) 
2576:  
2577: !      DEALLOCATE (R, QTRN, V, W, F, T, P, T4, PHI, FSITE, MST, KEP, PEP, SUMEPT, STAT = STATUS) 
2578:  
2579:       call RUN_INFO() 
2580:  
2581: !      END PROGRAM LWOTP_MD_NVE 
2582:       END SUBROUTINE RIGIDMD_NVE_CHANGE_DENSITY 
2583:  
2584: !     ------------------------------------------------------------------------------------ 
2585:   
2586:       SUBROUTINE NVECHDSTEP1 (S, PS, FCTMQT ) 
2587:  
2588:       IMPLICIT NONE 
2589:  
2590:       DOUBLE PRECISION, INTENT(IN) :: S, PS 
2591:       DOUBLE PRECISION, INTENT(IN) :: FCTMQT 
2592:  
2593:       END SUBROUTINE NVECHDSTEP1 
2594:  
2595: !     ------------------------------------------------------------------------------------ 
2596:   
2597:       SUBROUTINE NVECHDSTEP2 (PE) 
2598:  
2599:       !STEP 2 in the MOVE algorithm 
2600:       ! STEP 2: update V, P, PS 
2601:       ! PS = PS - PE*HALFDT   !! PE=potential energy 
2602:       ! loop over molecules 
2603:       !   update linear momentum (or velocity) 
2604:       !   V(I,:) = V(I,:) + S * F(I,:) * HALFDT / MASS    
2605:       !   update rotational momentum  
2606:       !   P(I,:) = P(I,:) +  2*S * DELT/2 * MATMUL(MX,T4(I,:)) 
2607:       ! end loop over molecules 
2608:        
2609:       USE MDCOMMONS, only : NMOL, HALFDT, MASS, QTRN, DELT, T4, P, V, F 
2610:  
2611:       IMPLICIT NONE 
2612:  
2613:       INTEGER :: I 
2614:       DOUBLE PRECISION :: MX(4,4), Q(4) 
2615:       DOUBLE PRECISION, INTENT(IN) :: PE 
2616:  
2617:       !PS     = PS - PE * HALFDT 
2618:  
2619:       DO I = 1, NMOL 
2620:  
2621:          V(I,:) = V(I,:) + F(I,:) * HALFDT / MASS           
2622:  
2623:          Q      = QTRN(I,:) 
2624:  
2625:          CALL CONSTRUCT_MX( MX, Q ) 
2626:  
2627:          P(I,:)     = P(I,:) +  DELT * MATMUL(MX,T4(I,:)) 
2628:  
2629:       ENDDO 
2630:  
2631:       END SUBROUTINE NVECHDSTEP2 
2632:  
2633: !     ------------------------------------------------------------------------------------ 
2634:   
2635:       SUBROUTINE NVECHDSTEP3 () 
2636:  
2637:       !STEP 3 in the MOVE algorithm 
2638:       ! STEP 3: update QTRN, P, PS with contributions from principle axis 3 
2639:       ! ZETASUM = 0 
2640:       ! loop over molecules 
2641:       !   ZETA(I,3) = dot_product( P(I,:) , MX(:,4) ) / ( 4 * S * JMI(3) ) 
2642:       !   ZETASUM = ZETASUM + ZETA(I,3)**2 
2643:       !   Q(I) = cos( ZETA(I,3)*HALFDT )*Q(I) + sin( ZETA(I,3)*HALFDT )*MX(:,4) 
2644:       !   P(I) = cos( ZETA(I,3)*HALFDT )*P(I) + sin( ZETA(I,3)*HALFDT )*MP(:,4) 
2645:       !   !   in the above, MP(:,:) is the analogous matrix to MX for vector P(I,:), 
2646:       !   !   i.e. MP(:,4) = (/-P(I,4),P(I,3),-P(I,2),P(I,1)/) 
2647:       ! end loop over molecules 
2648:       ! PS = PS + ZETASUM*JMI(3)*DELT 
2649:       ! contribution from principle axis 2 
2650:        
2651:       USE MDCOMMONS, only : QTRN, P, HALFDT, JMI, NMOL  
2652:  
2653:       IMPLICIT NONE 
2654:  
2655:       INTEGER :: I 
2656:       DOUBLE PRECISION :: SMPHSQ, Q(4), DQ(4), PHIT 
2657:       DOUBLE PRECISION :: PI(4), CP, SP 
2658:  
2659:       SMPHSQ = 0.D0 
2660:       DO I = 1, NMOL 
2661:  
2662:          Q         = QTRN(I,:) 
2663:          PI        = P(I,:) 
2664:          !DQ        = MX(:,4) 
2665:          DQ        = (/-Q(4),Q(3),-Q(2),Q(1)/) 
2666:          PHIT      = 0.25D0*DOT_PRODUCT(PI,DQ)*HALFDT/(JMI(3)) 
2667:          !SMPHSQ    = SMPHSQ + PHIT * PHIT 
2668:          CP        = COS(PHIT) 
2669:          SP        = SIN(PHIT) 
2670:          P(I,:)    = CP*PI + SP*(/-PI(4),PI(3),-PI(2),PI(1)/) 
2671:          QTRN(I,:) = CP*Q + SP*DQ 
2672:  
2673:       ENDDO 
2674:  
2675:       !PS     = PS + JMI(3)*SMPHSQ*DELT 
2676:  
2677:       END SUBROUTINE NVECHDSTEP3 
2678:  
2679: !     ------------------------------------------------------------------------------------ 
2680:   
2681:       SUBROUTINE NVECHDSTEP4 () 
2682:  
2683:       !STEP 4 in the MOVE algorithm 
2684:       ! STEP 4: update QTRN, P, PS with contributions from principle axis 2 
2685:       ! ZETASUM = 0 
2686:       ! loop over molecules 
2687:       !   ZETA(I,2) = dot_product( P(I,:) , MX(:,4) ) / ( 4 * S * JMI(2) ) 
2688:       !   ZETASUM = ZETASUM + ZETA(I,2)**2 
2689:       !   Q(I) = cos( ZETA(I,2)*HALFDT )*Q(I) + sin( ZETA(I,2)*HALFDT )*MX(:,2) 
2690:       !   P(I) = cos( ZETA(I,2)*HALFDT )*P(I) + sin( ZETA(I,2)*HALFDT )*MP(:,2) 
2691:       !   !   in the above, MP(:,:) is the analogous matrix to MX for vector P(I,:), 
2692:       !   !   i.e. MP(:,2) = (/-P(I,3),-P(I,4),P(I,1),P(I,2)/) 
2693:       ! end loop over molecules 
2694:       ! PS = PS + ZETASUM*JMI(2)*DELT 
2695:        
2696:       USE MDCOMMONS, only : QTRN, P, HALFDT, JMI, NMOL  
2697:  
2698:       IMPLICIT NONE 
2699:  
2700:       INTEGER :: I 
2701:       DOUBLE PRECISION :: SMPHSQ, Q(4), DQ(4), PHIT 
2702:       DOUBLE PRECISION :: PI(4), CP, SP 
2703:  
2704:       SMPHSQ = 0.D0 
2705:       DO I = 1, NMOL 
2706:  
2707:          Q         = QTRN(I,:) 
2708:          PI        = P(I,:) 
2709:          DQ        = (/-Q(3),-Q(4),Q(1),Q(2)/) 
2710:          PHIT      = 0.25D0*DOT_PRODUCT(PI,DQ)*HALFDT/(JMI(2)) 
2711:          SMPHSQ    = SMPHSQ + PHIT*PHIT 
2712:          CP        = COS(PHIT) 
2713:          SP        = SIN(PHIT) 
2714:          P(I,:)    = CP*PI + SP*(/-PI(3),-PI(4),PI(1),PI(2)/) 
2715:          QTRN(I,:) = CP*Q + SP*DQ 
2716:  
2717:       ENDDO 
2718:  
2719:       !PS     = PS + JMI(2)*SMPHSQ*DELT 
2720:  
2721:       END SUBROUTINE NVECHDSTEP4 
2722:  
2723: !     ------------------------------------------------------------------------------------ 
2724:   
2725:       SUBROUTINE NVECHDSTEP5 () 
2726:  
2727:       !STEP 5 in the MOVE algorithm 
2728:       ! STEP 5: update QTRN, P, PS with contributions from principle axis 1 
2729:       ! STEP 5 is analogous to STEP 3 and STEP 4 but for ZETA(:,1), i.e. the 
2730:       ! first principle axis, 
2731:       ! EXCEPT, the update of PS is different and involves the translational 
2732:       ! kinetic energy, the temperature (TMPFIX), the number of degrees of 
2733:       ! freedom (g), and HNOT, 
2734:        
2735:       USE MDCOMMONS, only : QTRN, P, HALFDT, JMI, NMOL, DELT, R, V 
2736:  
2737:       IMPLICIT NONE 
2738:  
2739:       INTEGER :: I 
2740:       DOUBLE PRECISION :: SMPHSQ, Q(4), DQ(4), PHIT 
2741:       DOUBLE PRECISION :: PI(4), CP, SP 
2742:       !DOUBLE PRECISION :: TKE !translational kinetic energy 
2743:  
2744:       SMPHSQ = 0.D0 
2745:       !TKE = 0 
2746:       DO I = 1, NMOL 
2747:  
2748:          R(I,:)    = R(I,:) + V(I,:) * DELT 
2749:          Q         = QTRN(I,:) 
2750:          PI        = P(I,:) 
2751:          DQ        = (/-Q(2),Q(1),Q(4),-Q(3)/) 
2752:          PHIT      = 0.25D0*DOT_PRODUCT(PI,DQ)*DELT/(JMI(1)) 
2753:          !SMPHSQ    = SMPHSQ + PHIT * PHIT 
2754:          CP        = COS(PHIT) 
2755:          SP        = SIN(PHIT) 
2756:          P(I,:)    = CP*PI + SP*(/-PI(2),PI(1),PI(4),-PI(3)/) 
2757:          QTRN(I,:) = CP*Q  + SP*DQ 
2758:          !TKE       = TKE + V(I,1)*V(I,1) + V(I,2)*V(I,2) + V(I,3)*V(I,3) 
2759:  
2760:       ENDDO 
2761:  
2762:       !PS = PS + (0.5D0*MASS*TKE/(S*S) + 2.D0*JMI(1)*SMPHSQ - DBLE(G)*TMPFIX*LOG(S) & 
2763:          !+ HNOT - DBLE(G)*TMPFIX) * DELT 
2764:  
2765:       END SUBROUTINE NVECHDSTEP5 
2766:  
2767: !     ---------------------------------------------------------------------------------------------- 
2768:  
2769:       SUBROUTINE NVECHDSCLVEL(TMPT, TMPR) 
2770:  
2771: ! scale linear velocities??? 
2772:  
2773:       USE MDCOMMONS, only : TMPFIX, V, W, QTRN, NMOL, JMI, P 
2774:  
2775:       IMPLICIT NONE 
2776:  
2777:       INTEGER          :: I 
2778:       DOUBLE PRECISION :: SCLTR, SCLRT, TMPT, TMPR 
2779:       DOUBLE PRECISION :: Q(4), RMXT(3,3), MX(4,4), W4(4) 
2780:  
2781:       SCLTR   = DSQRT(TMPFIX / TMPT) 
2782:       SCLRT   = DSQRT(TMPFIX / TMPR) 
2783:  
2784:       DO I = 1, NMOL 
2785:  
2786:          V(I,:) = V(I,:) * SCLTR 
2787:  
2788:          W(I,:) = W(I,:) * SCLRT 
2789:  
2790:          Q    = QTRN(I,:) 
2791:          CALL CONSTRUCT_RMXT(RMXT,Q) 
2792:  
2793: !     CONSTRUCT THE MATRIX MX FROM THE QUATERNION Q 
2794:  
2795:          CALL CONSTRUCT_MX(MX,Q) 
2796:  
2797:          W4(1)  = 0.D0 
2798:          W4(2)  = 2.D0*JMI(1)*W(I,1) 
2799:          W4(3)  = 2.D0*JMI(2)*W(I,2) 
2800:          W4(4)  = 2.D0*JMI(3)*W(I,3) 
2801:  
2802:          P(I,:) = MATMUL(MX,W4) 
2803:  
2804:       ENDDO 
2805:  
2806:       END SUBROUTINE NVECHDSCLVEL 
2807: ! 
2808: ! svn info  
2809: ! $Date: 2011-10-16 20:53:28 +0100 (Sun, 16 Oct 2011) $ 
2810: ! $Rev: 61 $ 
2811: ! 
2812:  
2813: !     ---------------------------------------------------------------------------------------------- 
2814:  
2815: !      PROGRAM LWOTP_MD_NVT 
2816:       SUBROUTINE RIGIDMD_NVT 
2817:       USE MDCOMMONS 
2818:  
2819:       IMPLICIT NONE 
2820:  
2821: ! IDUMP1: dump files.  always:  
2822: !                     enrg1.dat: KEPP, PEPP, EPP: kinetic, potential, energy per molecule 
2823: !                     tmpprs1.dat: temp, pressure? 
2824: !                     avenrg1.dat: average KE, PE, E 
2825: !                     avtmpprs1.dat: average temp, pressure? 
2826: !                     avtmptr1.dat: average tranlational, rotational temperature 
2827: ! IDUMP2: dump files.  ISTEP > NEQ 
2828: !                     pos1.dat:  xyz positions of CoM of molecules 
2829: !                     ortn1.dat: quaternions giving the orientation of molecules 
2830: !                     rot1.dat:  phi1, phi2, phi3 ???? 
2831: !                     efm1.dat:  energy fluctuation metric: time from NEQ, and variance over molecules of time averaged energy  
2832: ! IDUMP3: CALL SCLVEL.  ISTEP <= NBATH: scale linear velocities ??? 
2833: ! IDUMP4: dump files.  ISTEP > NEQ:  
2834: !                     blockav1.dat: average PE, temp, pressure 
2835: !                     lvel1.dat:    linear velocities 
2836: !                     avel1.dat:    angular velocities 
2837: ! TKE: translational kinetic energy 
2838: ! RKE: rotational kinetic energy 
2839: ! MX:  is a representation of the quaternion as a matrix in such a way 
2840: !      that quaternion addition and multiplication correspond to matrix 
2841: !      addition and matrix multiplication 
2842: !      In this way the conjugate of the quaternion corresponds to the 
2843: !      transpose of the matrix.  The fourth power of the norm of a 
2844: !      quaternion is the determinant of the corresponding matrix. 
2845: !      Complex numbers are block diagonal matrices with two 2x2 blocks. 
2846: ! EFMT: Energy Fluctuation Metric: the variance over molecules of the 
2847: !       time averaged energy 
2848: ! VRL: accumulates the forces dotted into the positions Fij.dot(Ri - Rj) which 
2849: !      is used to calculate the pressure 
2850:       INTEGER          :: NSTEP, NBATH, NEQ, ISTEP, ICNT, IDUMP1, IDUMP2, IDUMP4, I!, K!, STATUS 
2851:       DOUBLE PRECISION :: RCUT6, RCUT12, TMP, TMPTR, TMPRT, RHO, VLM 
2852:       DOUBLE PRECISION :: Q(4), RMXT(3,3)!, MX(4,4), W4(4) 
2853:       DOUBLE PRECISION :: PRS, SUMPRS, AVPRS, PE, KE, TKE, RKE, PEPP, KEPP, EPP, VRL 
2854:       DOUBLE PRECISION :: SUME, SUMKE, SUMPE, SUMTMT, SUMTMR, SUMTMP 
2855:       DOUBLE PRECISION :: AVE, AVKE, AVPE, AVTMPT, AVTMPR, AVTMP, EFMT, TIME 
2856:       DOUBLE PRECISION :: BSUMPE, BSUMTMP, BSUMPRS, AVPEB, AVTMPB, AVPRSB, DBLE 
2857:       DOUBLE PRECISION :: MQ=0.1D0, DH=1.D0 !only NVT 
2858:       DOUBLE PRECISION :: TMPINT !initial temperature 
2859:       logical :: bool1, bool2, loadHNOT 
2860:       !double precision :: rand 
2861:       !CALL RANDOM_SEED() 
2862:       !CALL RANDOM_number( rand ) 
2863:       !write(*,*) rand 
2864:       !DOUBLE PRECISION :: dummyvec3(3), dummyvec4(4) 
2865:  
2866:  
2867:       OPEN (UNIT = 1, FILE = 'parameternvt1.inp', STATUS = 'UNKNOWN') 
2868:  
2869: !     READ INPUT PARAMETERS 
2870:  
2871:       READ (1, *) 
2872:       READ (1, *) NMOL, NSITE 
2873:       READ (1, *)  
2874:       READ (1, *) RCUT, RLIST 
2875:       READ (1, *) 
2876:       READ (1, *) TMPFIX 
2877:       READ (1, *) 
2878:       READ (1, *) RHO 
2879:       READ (1, *) 
2880:       READ (1, *) DELT 
2881:       READ (1, *) 
2882:       READ (1, *) S, PS, MQ 
2883:       READ (1, *)  
2884:       READ (1, *) NSTEP, NBATH, NEQ 
2885:       READ (1, *)  
2886:       READ (1, *) IDUMP1, IDUMP2, IDUMP4 
2887:   
2888:       CLOSE (1) 
2889:       TMPINT = TMPFIX !initial temperature 
2890:  
2891: !     CALCULATE FROM INPUT PARAMETERS 
2892:  
2893:       G      = 6 * NMOL - 3 
2894:       NTST   = NMOL * NSITE             ! TOTAL NUMBER OF SITES 
2895:       RCUTSQ = RCUT * RCUT 
2896:       RLSTSQ = RLIST * RLIST 
2897:       EPS4   = 4.D0  ! 4*epsilon 
2898:       RCUT6  = (1.D0 / RCUT) ** 6 
2899:       RCUT12 = RCUT6 * RCUT6 
2900:       CNSTA  = (6.D0 * RCUT12 - 3.D0 * RCUT6) / RCUTSQ 
2901:       CNSTB  = 4.D0 * RCUT6 - 7.D0 * RCUT12 
2902:       HALFDT = 0.5D0 * DELT 
2903:       FCTMQT = 0.5D0*HALFDT/MQ !used only in NVT 
2904:   
2905: !     ALLOCATE THE ARRAYS 
2906:  
2907: !      ALLOCATE (R(NMOL,3), QTRN(NMOL,4), V(NMOL,3), W(NMOL,3), F(NMOL,3), T(NMOL,3), P(NMOL,4), & 
2908: !      T4(NMOL,4), PHI(NMOL,3), FSITE(NTST,3), MST(NSITE,3), KEP(NMOL), PEP(NMOL), SUMEPT(NMOL), & 
2909: !      STAT = STATUS) 
2910:  
2911: !      IF (STATUS /= 0) STOP 'ALLOCATION PROBLEM' 
2912:   
2913: !     BOX LENGTH 
2914:  
2915:       VLM    = DBLE(NMOL) / RHO 
2916:       BOXL   = VLM ** (1.D0/3.D0) 
2917:  
2918:       CALL DEFMOL() 
2919:  
2920: !      WRITE(*,*) RCUT, RCUT2, RLIST 
2921:       WRITE(*,*) "Read in input parameters" 
2922:  
2923:       INQUIRE(FILE="initpos1.dat", EXIST=bool1 )  
2924:       INQUIRE(FILE="initortn1.dat", EXIST=bool2 )  
2925:       if ( bool1 .and. bool2 ) then 
2926:         !start from an existing configuration 
2927:         !CALL EQCON() 
2928:         write(*,*) "reading init files: initpos1.dat and initortn1.dat" 
2929:         CALL LOAD_POS_ORTN() 
2930:       else 
2931:         !     start from scratch 
2932:         !     START WITH A SIMPLE CUBIC LATTICE CONFIGURATION 
2933:         write(*,*) "starting from a simple cubic lattice" 
2934:         CALL RIGIDMD_SC() 
2935:         !     SET INITIAL ORIENTATIONS 
2936:         CALL INTORN() 
2937:       endif 
2938:  
2939: !     SET INITIAL LINEAR and angular VELICITIES 
2940:       INQUIRE(FILE="initlinvel1.dat", EXIST=bool1 )  
2941:       INQUIRE(FILE="initangvel1.dat", EXIST=bool2 )  
2942:       if ( bool1 .and. bool2 ) then 
2943:         write(*,*) "reading init files: initlinvel1.dat and initangvel1.dat" 
2944:         CALL LOAD_LIN_ANG_VEL() 
2945:       else 
2946:         write(*,*) "assigning random velocities" 
2947:         CALL INTVEL(TMPINT) 
2948:       endif 
2949:       V=V*S !convert from real to virtual velocity 
2950:  
2951: !     SET INITIAL THERMOSTAT VALUES 
2952:       INQUIRE(FILE="initthermostat1.dat", EXIST=bool1 )  
2953:       if ( bool1 ) then 
2954:         write(*,*) "reading init file: initthermostat1.dat" 
2955:         OPEN (UNIT=17, FILE='initthermostat1.dat', STATUS='UNKNOWN') 
2956:         read(17,*) S, PS, HNOT 
2957:         CLOSE(17) 
2958:         loadHNOT = .true. 
2959:       ELSE 
2960:         !S = 1.D0 
2961:         !PS = 0.D0 
2962:         write(*,*) "setting thermostat initial values S=",S, "PS=",PS 
2963:         loadHNOT = .false. 
2964:       endif 
2965:  
2966:  
2967: !     START SIMULATION 
2968:  
2969:       ISTEP  = 0   
2970:       ICNT   = 0 
2971:       SUME   = 0.D0 
2972:       SUMKE  = 0.D0 
2973:       SUMPE  = 0.D0 
2974:       SUMTMT = 0.D0 
2975:       SUMTMR = 0.D0 
2976:       SUMTMP = 0.D0 
2977:       SUMPRS = 0.D0 
2978:        
2979:       PHI(:,:) = 0.D0 
2980:       KE       = 0.D0 
2981:  
2982:  
2983:       UPDATE = .TRUE. 
2984:  
2985: !     CALCULATE FORCES AND TORQUES AT TIME t=0 
2986:       CALL FORQUE(PE, VRL) 
2987:  
2988:       DO I = 1, NMOL 
2989:  
2990: !        TORQUE AT TIME t=0 IN THE BODY FRAME USING THE TRANSPOSE OF THE ROTATION MATRIX 
2991:          Q    = QTRN(I,:) 
2992:          CALL CONSTRUCT_RMXT( RMXT, Q ) 
2993:          T4(I,1)   = 0.D0 
2994:          T4(I,2:4) = MATMUL(RMXT,T(I,:)) 
2995:          ! intitialize P, the angular momentum.  P is virtual, W is real 
2996:          !CALL CONSTRUCT_MX( MX, Q ) 
2997:          !W4(1)  = 0.D0 
2998:          !W4(2)  = 2.D0*W(I,1)*JMI(1) 
2999:          !W4(3)  = 2.D0*W(I,2)*JMI(2) 
3000:          !W4(4)  = 2.D0*W(I,3)*JMI(3) 
3001:          !P(I,:) = MATMUL(MX,W4*S) 
3002:          call PfromW( P(I,:), W(I,:), JMI, Q, S ) 
3003:          !if ( I .eq. 1 ) then 
3004:            !dummyvec4=P(I,:) 
3005:            !dummyvec3=W(I,:) 
3006:            !write(*,*) dummyvec4 
3007:            !write(*,*) dummyvec3 
3008:            !call WfromP( dummyvec4, dummyvec3, JMI, Q, S ) 
3009:            !write(*,*)   
3010:            !write(*,*) dummyvec4 
3011:            !write(*,*) dummyvec3 
3012:          !endif 
3013:   
3014:          !calculate the energy 
3015:          KE = KE + 0.5D0*MASS*(V(I,1)*V(I,1) + V(I,2)*V(I,2) + V(I,3)*V(I,3))/S**2 & 
3016:                  + 0.5D0*(W(I,1)*W(I,1)*JMI(1) + W(I,2)*W(I,2)*JMI(2) + W(I,3)*W(I,3)*JMI(3)) ! here V is virtual and W is real 
3017:  
3018:       ENDDO 
3019:  
3020:       IF ( loadHNOT .eqv. .false. ) THEN 
3021:         HNOT = KE + PE + 0.5D0*PS*PS/MQ + DBLE(G)*TMPFIX*LOG(S) 
3022:       ENDIF 
3023:  
3024:       !print initial values for bug testing only 
3025:       IF ( .true. ) then 
3026:         OPEN (UNIT=71, FILE='pos0.dat', STATUS='UNKNOWN') 
3027:         OPEN (UNIT=72, FILE='ortn0.dat', STATUS='UNKNOWN') 
3028:         OPEN (UNIT=73, FILE='lvel0.dat', STATUS='UNKNOWN') 
3029:         OPEN (UNIT=74, FILE='avel0.dat', STATUS='UNKNOWN') 
3030:         WRITE(*,*) "Writing to initial output files" 
3031:         do I=1,NMOL 
3032:           WRITE (71, *)  R(I,1), R(I,2), R(I,3) 
3033:           WRITE (72, 900)  QTRN(I,1), QTRN(I,2), QTRN(I,3), QTRN(I,4)  
3034:           WRITE (73, *)  V(I,1)/S, V(I,2)/S, V(I,3)/S 
3035:           WRITE (74, *) W(I,1), W(I,2), W(I,3) 
3036:         enddo 
3037:         close(71) 
3038:         close(72) 
3039:         close(73) 
3040:         close(74) 
3041:         OPEN (UNIT = 75, FILE = 'thermostat0.dat', STATUS = 'UNKNOWN') 
3042:           WRITE (75, *) S, PS, HNOT 
3043:         CLOSE(UNIT=75) 
3044:            !OPEN (UNIT=29, FILE='angmom0.dat', STATUS='UNKNOWN', POSITION='APPEND') 
3045:            !DO I = 1, NMOL 
3046:              !WRITE (29, 900)  P(I,1), P(I,2), P(I,3), P(I,4)  
3047:            !ENDDO 
3048:            !CLOSE (UNIT=29,  STATUS='KEEP')   
3049:       endif 
3050:  
3051:  
3052:       DO WHILE (ISTEP < NSTEP)             
3053:  
3054:          ISTEP = ISTEP + 1 
3055:          ICNT  = ICNT + 1 
3056:          IF(MOD(ISTEP,1000).EQ.0) WRITE(*,*) "Step number:", ISTEP, "of ", NSTEP 
3057: !     ADVANCE POSITIONS, ORIENTATIONS, AND THEIR TIME DERIVATIVES 
3058: !     MOVE also calls CHECK() and FORQUE 
3059:          
3060:          CALL NVTMOVE(ISTEP, NEQ, PE, VRL, TKE, RKE) 
3061:  
3062: !     CALCULATE TEMPERATURES, PRESSURE AND TOTEL ENERGY AT TIME t 
3063:  
3064:          TMPTR   = 2.D0 * TKE  / DBLE(3 * NMOL - 3) 
3065:          TMPRT   = 2.D0 * RKE / DBLE(3 * NMOL)            
3066:          KE      = TKE + RKE 
3067:          !DH is used only in NVT 
3068:          DH      = ((KE + PE + 0.5D0*PS*PS/MQ + DBLE(G)*TMPFIX*LOG(S) - HNOT)/HNOT)  
3069:          TMP     = 2.D0 * KE / DBLE(G) 
3070:          PEPP    = PE / DBLE(NMOL) 
3071:          KEPP    = KE / DBLE(NMOL) 
3072:          EPP     = PEPP + KEPP 
3073:          PRS     = (2.D0 * TKE + VRL) / (3.D0 * VLM) 
3074:          SUME    = SUME + EPP 
3075:          SUMKE   = SUMKE + KEPP 
3076:          SUMPE   = SUMPE + PEPP 
3077:          SUMTMT  = SUMTMT + TMPTR 
3078:          SUMTMR  = SUMTMR + TMPRT 
3079:          SUMTMP  = SUMTMP + TMP 
3080:          SUMPRS  = SUMPRS + PRS 
3081:  
3082:          !do block averaging of potential energy, temperature, and pressure 
3083:          IF (ISTEP > NEQ) THEN 
3084:  
3085:             BSUMPE  = BSUMPE + PEPP 
3086:             BSUMTMP = BSUMTMP + TMP 
3087:             BSUMPRS = BSUMPRS + PRS 
3088:  
3089:             IF (MOD((ISTEP-NEQ),IDUMP4) == 0) THEN 
3090:  
3091:                AVPEB  = BSUMPE / DBLE(IDUMP4) 
3092:                AVTMPB = BSUMTMP / DBLE(IDUMP4) 
3093:                AVPRSB = BSUMPRS / DBLE(IDUMP4)   
3094:  
3095:                OPEN (UNIT=25, FILE='blockav1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
3096:                  WRITE (25, *) AVPEB, AVTMPB, AVPRSB 
3097:                CLOSE (UNIT=25, STATUS='KEEP') 
3098:  
3099:                BSUMPE  = 0.D0 
3100:                BSUMTMP = 0.D0 
3101:                BSUMPRS = 0.D0 
3102:  
3103:             ENDIF 
3104:  
3105:          ENDIF 
3106:   
3107:          !write out energies, temperatures, and their averages, etc. 
3108:          IF (MOD(ISTEP, IDUMP1) .EQ. 0) THEN 
3109:  
3110:             AVE    = SUME / DBLE(ICNT) 
3111:             AVKE   = SUMKE / DBLE(ICNT) 
3112:             AVPE   = SUMPE / DBLE(ICNT) 
3113:             AVTMPT = SUMTMT / DBLE(ICNT) 
3114:             AVTMPR = SUMTMR / DBLE(ICNT) 
3115:             AVTMP  = SUMTMP / DBLE(ICNT) 
3116:             AVPRS  = SUMPRS / DBLE(ICNT) 
3117:  
3118:             OPEN (UNIT=3, FILE='enrg1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
3119:                  WRITE (3, *) KEPP, PEPP, EPP  
3120:             CLOSE (UNIT=3, STATUS='KEEP') 
3121:             OPEN (UNIT=4, FILE='tmpprs1.dat', STATUS = 'UNKNOWN', POSITION='APPEND') 
3122:                  WRITE (4, *) TMP, PRS 
3123:             CLOSE (UNIT=4, STATUS='KEEP')  
3124:  
3125:             OPEN (UNIT=33, FILE='avenrg1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
3126:                  WRITE (33, *) AVKE, AVPE, AVE 
3127:             CLOSE (UNIT=33, STATUS='KEEP') 
3128:             OPEN (UNIT=34, FILE='avtmpprs1.dat', STATUS = 'UNKNOWN', POSITION='APPEND') 
3129:                  WRITE (34, *) AVTMP, AVPRS 
3130:             CLOSE (UNIT=34, STATUS='KEEP')  
3131:  
3132:             OPEN (UNIT=35, FILE='avtmptr1.dat', STATUS = 'UNKNOWN', POSITION='APPEND') 
3133:                  WRITE (35, *) AVTMPT, AVTMPR 
3134:             CLOSE (UNIT=35, STATUS='KEEP') 
3135:  
3136:             OPEN (UNIT=36, FILE='error1.dat', STATUS = 'UNKNOWN', POSITION='APPEND') 
3137:                  WRITE (36, *) DH,  DH*HNOT+HNOT 
3138:             CLOSE (UNIT=36, STATUS='KEEP') 
3139:  
3140:             OPEN (UNIT=37, FILE='thermostat1.dat', STATUS = 'UNKNOWN', POSITION='APPEND') 
3141:                  WRITE (37, *) S, PS 
3142:             CLOSE (UNIT=37, STATUS='KEEP') 
3143:  
3144:          ENDIF 
3145:  
3146:          ! calculate energy averaged over time and molecules. 
3147:          ! calculate variance over molecules of time averaged energy 
3148:          IF (ISTEP > NEQ) THEN 
3149:  
3150:             CALL EFM (EFMT, TIME) 
3151:  
3152:          ENDIF 
3153:  
3154:          !write out positions, orientations, and angular change since t=NEQ 
3155:          IF (ISTEP > NEQ .AND. MOD(ISTEP, IDUMP2) == 0) THEN    
3156:  
3157:            OPEN (UNIT=7,  FILE='pos1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
3158:            OPEN (UNIT=8,  FILE='ortn1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
3159:            OPEN (UNIT=28, FILE='rot1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
3160:            OPEN (UNIT=29, FILE='efm1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
3161:  
3162:            DO I = 1, NMOL 
3163:               
3164:              WRITE (7, *)  R(I,1), R(I,2), R(I,3) 
3165:              WRITE (8, 900)  QTRN(I,1), QTRN(I,2), QTRN(I,3), QTRN(I,4)  
3166:              WRITE (28, *) PHI(I,1), PHI(I,2), PHI(I,3)   
3167:              900 FORMAT (1x,4ES25.16) 
3168:  
3169:            ENDDO 
3170:  
3171:            WRITE(29, *) TIME, EFMT 
3172:  
3173:  
3174:            CLOSE (UNIT=7,  STATUS='KEEP')   
3175:            CLOSE (UNIT=8,  STATUS='KEEP')   
3176:            CLOSE (UNIT=28, STATUS='KEEP') 
3177:            CLOSE (UNIT=29, STATUS='KEEP') 
3178:  
3179:            !OPEN (UNIT=29, FILE='angmom1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
3180:            !DO I = 1, NMOL 
3181:              !WRITE (29, 900)  P(I,1), P(I,2), P(I,3), P(I,4)  
3182:            !ENDDO 
3183:            !CLOSE (UNIT=29,  STATUS='KEEP')   
3184:  
3185:            IF (MOD(ISTEP, IDUMP4) .EQ. 0) THEN    
3186:  
3187:              OPEN (UNIT=9, FILE='lvel1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
3188:              OPEN (UNIT=10, FILE='avel1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
3189:  
3190:              DO I = 1, NMOL 
3191:  
3192:                WRITE (9, *)  V(I,1)/S, V(I,2)/S, V(I,3)/S 
3193:                WRITE (10, *) W(I,1), W(I,2), W(I,3) 
3194:  
3195:              ENDDO 
3196:  
3197:              CLOSE (UNIT = 9, STATUS = 'KEEP')   
3198:              CLOSE (UNIT = 10, STATUS ='KEEP')   
3199:  
3200:            ENDIF 
3201:  
3202:          ENDIF 
3203:                 
3204:          !reinitialize energy averages 
3205:          IF ( ISTEP == NEQ) THEN 
3206:                
3207:             ICNT   = 0 
3208:             SUME   = 0.D0 
3209:             SUMKE  = 0.D0 
3210:             SUMPE  = 0.D0 
3211:             SUMTMT = 0.D0 
3212:             SUMTMR = 0.D0 
3213:             SUMTMP = 0.D0 
3214:             SUMPRS = 0.D0         
3215:             TIME   = 0.D0 
3216:  
3217:             BSUMPE  = 0.D0 
3218:             BSUMTMP = 0.D0 
3219:             BSUMPRS = 0.D0 
3220:                           
3221:          ENDIF   
3222:  
3223:       ENDDO 
3224:  
3225:       !program ended, print out the final data 
3226:  
3227:       OPEN (UNIT = 21, FILE = 'finalpos1.dat', STATUS = 'UNKNOWN') 
3228:       OPEN (UNIT = 22, FILE = 'finalortn1.dat', STATUS = 'UNKNOWN') 
3229:       OPEN (UNIT = 23, FILE = 'finallinvel1.dat', STATUS = 'UNKNOWN') 
3230:       OPEN (UNIT = 24, FILE = 'finalangvel1.dat', STATUS = 'UNKNOWN') 
3231:        
3232:       DO I = 1, NMOL 
3233:  
3234:          WRITE (21, *) R(I,1), R(I,2), R(I,3) 
3235:          WRITE (22, 900) QTRN(I,1), QTRN(I,2), QTRN(I,3), QTRN(I,4) 
3236:          WRITE (23, *) V(I,1)/S, V(I,2)/S, V(I,3)/S 
3237:          WRITE (24, *) W(I,1), W(I,2), W(I,3) 
3238:  
3239:       ENDDO 
3240:  
3241:       CLOSE(UNIT=21) 
3242:       CLOSE(UNIT=22) 
3243:       CLOSE(UNIT=23) 
3244:       CLOSE(UNIT=24) 
3245:  
3246:       OPEN (UNIT = 25, FILE = 'finalthermostat1.dat', STATUS = 'UNKNOWN') 
3247:          WRITE (25, *) S, PS, HNOT 
3248:       CLOSE(UNIT=25) 
3249:            !OPEN (UNIT=29, FILE='finalangmom1.dat', STATUS='UNKNOWN', POSITION='APPEND') 
3250:            !DO I = 1, NMOL 
3251:              !WRITE (29, 900)  P(I,1), P(I,2), P(I,3), P(I,4)  
3252:            !ENDDO 
3253:            !CLOSE (UNIT=29,  STATUS='KEEP')   
3254:  
3255: !      DEALLOCATE (R, QTRN, V, W, F, T, P, T4, PHI, FSITE, MST, KEP, PEP, SUMEPT, STAT = STATUS) 
3256:  
3257:       call RUN_INFO() 
3258:  
3259: !      END PROGRAM LWOTP_MD_NVT 
3260:       END SUBROUTINE RIGIDMD_NVT 
3261: !     ------------------------------------------------------------------------------------ 
3262:   
3263:       SUBROUTINE NVTSTEP1 (S, PS, FCTMQT ) 
3264:  
3265:       !STEP 1 in the MOVE algorithm 
3266:       ! STEP 1: update S, PS 
3267:       ! S = S*(1.D0 + PS*0.5D0*HALFDT/MQ)**2 
3268:       ! PS = PS/(1.D0 + PS*0.5D0*HALFDT/MQ) 
3269:  
3270:       IMPLICIT NONE 
3271:  
3272:       DOUBLE PRECISION, INTENT(INOUT) :: S, PS 
3273:       DOUBLE PRECISION, INTENT(IN) :: FCTMQT 
3274:       DOUBLE PRECISION :: FCTR 
3275:  
3276:       FCTR = 1.D0 + FCTMQT*PS 
3277:       S    = S*FCTR*FCTR 
3278:       PS   = PS/FCTR     
3279:  
3280:       END SUBROUTINE NVTSTEP1 
3281:  
3282: !     ------------------------------------------------------------------------------------ 
3283:   
3284:       SUBROUTINE NVTSTEP2 (PE) 
3285:  
3286:       !STEP 2 in the MOVE algorithm 
3287:       ! STEP 2: update PS, V, and P from the Torque 
3288:       ! PS = PS - PE*HALFDT   !! PE=potential energy 
3289:       ! loop over molecules 
3290:       !   update linear momentum (or velocity) 
3291:       !   V(I,:) = V(I,:) + S * F(I,:) * HALFDT / MASS    
3292:       !   update rotational momentum  
3293:       !   P(I,:) = P(I,:) +  2*S * DELT/2 * MATMUL(MX,T4(I,:)) 
3294:       ! end loop over molecules 
3295:        
3296:       USE MDCOMMONS, only : PS, HALFDT, DELT, NMOL, V, QTRN, S, F, MASS, P, T4 
3297:  
3298:       IMPLICIT NONE 
3299:  
3300:       INTEGER :: I 
3301:       DOUBLE PRECISION :: MX(4,4), Q(4) 
3302:       DOUBLE PRECISION, INTENT(IN) :: PE 
3303:  
3304:       PS     = PS - PE * HALFDT 
3305:  
3306:       DO I = 1, NMOL 
3307:  
3308:          V(I,:) = V(I,:) + S * F(I,:) * HALFDT / MASS           
3309:  
3310:          Q      = QTRN(I,:) 
3311:  
3312:          CALL CONSTRUCT_MX( MX, Q ) 
3313:  
3314:          P(I,:)     = P(I,:) +  S * DELT * MATMUL(MX,T4(I,:)) 
3315:  
3316:       ENDDO 
3317:  
3318:       END SUBROUTINE NVTSTEP2 
3319:  
3320: !     ------------------------------------------------------------------------------------ 
3321:   
3322:       SUBROUTINE NVTSTEP3 () 
3323:  
3324:       !STEP 3 in the MOVE algorithm 
3325:       ! STEP 3: update QTRN, P, PS with contributions from principle axis 3 
3326:       ! ZETASUM = 0 
3327:       ! loop over molecules 
3328:       !   ZETA(I,3) = dot_product( P(I,:) , MX(:,4) ) / ( 4 * S * JMI(3) ) 
3329:       !   ZETASUM = ZETASUM + ZETA(I,3)**2 
3330:       !   Q(I) = cos( ZETA(I,3)*HALFDT )*Q(I) + sin( ZETA(I,3)*HALFDT )*MX(:,4) 
3331:       !   P(I) = cos( ZETA(I,3)*HALFDT )*P(I) + sin( ZETA(I,3)*HALFDT )*MP(:,4) 
3332:       !   !   in the above, MP(:,:) is the analogous matrix to MX for vector P(I,:), 
3333:       !   !   i.e. MP(:,4) = (/-P(I,4),P(I,3),-P(I,2),P(I,1)/) 
3334:       ! end loop over molecules 
3335:       ! PS = PS + ZETASUM*JMI(3)*DELT 
3336:       ! contribution from principle axis 2 
3337:        
3338:       USE MDCOMMONS, only : NMOL, QTRN, P, S, PS, JMI, HALFDT, DELT 
3339:  
3340:       IMPLICIT NONE 
3341:  
3342:       INTEGER :: I 
3343:       DOUBLE PRECISION :: SMPHSQ, Q(4), DQ(4), PHIT 
3344:       DOUBLE PRECISION :: PI(4), CP, SP 
3345:  
3346:       SMPHSQ = 0.D0 
3347:       DO I = 1, NMOL 
3348:  
3349:          Q         = QTRN(I,:) 
3350:          PI        = P(I,:) 
3351:          !DQ        = MX(:,4) 
3352:          DQ        = (/-Q(4),Q(3),-Q(2),Q(1)/) 
3353:          PHIT      = 0.25D0*DOT_PRODUCT(PI,DQ)/(S*JMI(3)) 
3354:          SMPHSQ    = SMPHSQ + PHIT * PHIT 
3355:          CP        = COS(PHIT*HALFDT) 
3356:          SP        = SIN(PHIT*HALFDT) 
3357:          P(I,:)    = CP*PI + SP*(/-PI(4),PI(3),-PI(2),PI(1)/) 
3358:          QTRN(I,:) = CP*Q + SP*DQ 
3359:  
3360:       ENDDO 
3361:  
3362:       PS     = PS + JMI(3)*SMPHSQ*DELT 
3363:  
3364:       END SUBROUTINE NVTSTEP3 
3365:  
3366: !     ------------------------------------------------------------------------------------ 
3367:   
3368:       SUBROUTINE NVTSTEP4 () 
3369:  
3370:       !STEP 4 in the MOVE algorithm 
3371:       ! STEP 4: update QTRN, P, PS with contributions from principle axis 2 
3372:       ! ZETASUM = 0 
3373:       ! loop over molecules 
3374:       !   ZETA(I,2) = dot_product( P(I,:) , MX(:,4) ) / ( 4 * S * JMI(2) ) 
3375:       !   ZETASUM = ZETASUM + ZETA(I,2)**2 
3376:       !   Q(I) = cos( ZETA(I,2)*HALFDT )*Q(I) + sin( ZETA(I,2)*HALFDT )*MX(:,2) 
3377:       !   P(I) = cos( ZETA(I,2)*HALFDT )*P(I) + sin( ZETA(I,2)*HALFDT )*MP(:,2) 
3378:       !   !   in the above, MP(:,:) is the analogous matrix to MX for vector P(I,:), 
3379:       !   !   i.e. MP(:,2) = (/-P(I,3),-P(I,4),P(I,1),P(I,2)/) 
3380:       ! end loop over molecules 
3381:       ! PS = PS + ZETASUM*JMI(2)*DELT 
3382:        
3383:       USE MDCOMMONS, only : NMOL, QTRN, P, S, PS, JMI, HALFDT, DELT 
3384:  
3385:       IMPLICIT NONE 
3386:  
3387:       INTEGER :: I 
3388:       DOUBLE PRECISION :: SMPHSQ, Q(4), DQ(4), PHIT 
3389:       DOUBLE PRECISION :: PI(4), CP, SP 
3390:  
3391:       SMPHSQ = 0.D0 
3392:       DO I = 1, NMOL 
3393:  
3394:          Q         = QTRN(I,:) 
3395:          PI        = P(I,:) 
3396:          DQ        = (/-Q(3),-Q(4),Q(1),Q(2)/) 
3397:          PHIT      = 0.25D0*DOT_PRODUCT(PI,DQ)/(S*JMI(2)) 
3398:          SMPHSQ    = SMPHSQ + PHIT*PHIT 
3399:          CP        = COS(PHIT*HALFDT) 
3400:          SP        = SIN(PHIT*HALFDT) 
3401:          P(I,:)    = CP*PI + SP*(/-PI(3),-PI(4),PI(1),PI(2)/) 
3402:          QTRN(I,:) = CP*Q + SP*DQ 
3403:  
3404:       ENDDO 
3405:  
3406:       PS     = PS + JMI(2)*SMPHSQ*DELT 
3407:  
3408:       END SUBROUTINE NVTSTEP4 
3409:  
3410: !     ------------------------------------------------------------------------------------ 
3411:   
3412:       SUBROUTINE NVTSTEP5 () 
3413:  
3414:       !STEP 5 in the MOVE algorithm 
3415:       ! STEP 5: update R QTRN, P, PS with contributions from principle axis 1 
3416:       ! STEP 5 is analogous to STEP 3 and STEP 4 but for ZETA(:,1), i.e. the 
3417:       ! first principle axis, 
3418:       ! EXCEPT, the update of PS is different and involves the translational 
3419:       ! kinetic energy, the temperature (TMPFIX), the number of degrees of 
3420:       ! freedom (g), and HNOT.  Also, R is updated 
3421:        
3422:       USE MDCOMMONS, only : NMOL, QTRN, P, S, PS, JMI, HALFDT, DELT, R, V, G, & 
3423:         TMPFIX, HNOT, MASS 
3424:  
3425:       IMPLICIT NONE 
3426:  
3427:       INTEGER :: I 
3428:       DOUBLE PRECISION :: SMPHSQ, Q(4), DQ(4), PHIT, DBLE 
3429:       DOUBLE PRECISION :: PI(4), CP, SP 
3430:       DOUBLE PRECISION :: TKE !translational kinetic energy 
3431:  
3432:       SMPHSQ = 0.D0 
3433:       TKE = 0 
3434:       DO I = 1, NMOL 
3435:  
3436:          R(I,:)    = R(I,:) + V(I,:) * DELT / S 
3437:          Q         = QTRN(I,:) 
3438:          PI        = P(I,:) 
3439:          DQ        = (/-Q(2),Q(1),Q(4),-Q(3)/) 
3440:          PHIT      = 0.25D0*DOT_PRODUCT(PI,DQ)/(S*JMI(1)) 
3441:          SMPHSQ    = SMPHSQ + PHIT * PHIT 
3442:          CP        = COS(PHIT*DELT) 
3443:          SP        = SIN(PHIT*DELT) 
3444:          P(I,:)    = CP*PI + SP*(/-PI(2),PI(1),PI(4),-PI(3)/) 
3445:          QTRN(I,:) = CP*Q  + SP*DQ 
3446:          TKE       = TKE + V(I,1)*V(I,1) + V(I,2)*V(I,2) + V(I,3)*V(I,3) 
3447:  
3448:       ENDDO 
3449:  
3450:       PS = PS + (0.5D0*MASS*TKE/(S*S) + 2.D0*JMI(1)*SMPHSQ & 
3451:            - DBLE(G)*TMPFIX*(1.D0+LOG(S)) + HNOT ) * DELT 
3452:  
3453:       END SUBROUTINE NVTSTEP5 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0