hdiff output

r30589/multipot.f90 2016-06-13 18:30:06.760981750 +0100 r30588/multipot.f90 2016-06-13 18:30:07.408990491 +0100
132: 132: 
133:         READ(22,*) DUMMYCHAR133:         READ(22,*) DUMMYCHAR
134:         READ(22,*) DUMMYCHAR134:         READ(22,*) DUMMYCHAR
135:         ! Now read in the atom numbers from the following line(s).135:         ! Now read in the atom numbers from the following line(s).
136:         ! The required input format will vary between potentials, although most will use the format described under CASE DEFAULT136:         ! The required input format will vary between potentials, although most will use the format described under CASE DEFAULT
137: 137: 
138:         SELECT CASE(POTTYPES(J1))138:         SELECT CASE(POTTYPES(J1))
139: 139: 
140:         ! For potentials which only operate between specific pairs of atoms, such as harmonic springs, each line should contain an interacting pair.140:         ! For potentials which only operate between specific pairs of atoms, such as harmonic springs, each line should contain an interacting pair.
141:         CASE('HSPR')141:         CASE('HSPR')
142:             DO J2 = 1, NATOM_BY_POT(J1)/2  ! J2 ranges over the number of springs (Note, NATOM_BY_POT contains the number of atoms, not the number of pairs)142:             DO J2 = 1, NATOM_BY_POT(J1)/2  ! J2 ranges over the number of springs (Note, NATOM_BY_POT contains the actual no of atoms)
143:                 ! We only want to compute the potential once for each spring.143:                 ! We only want to compute the potential once for each spring.
144: 144: 
145:                 READ(22,*) ATOM1, ATOM2145:                 READ(22,*) ATOM1, ATOM2
146: 146: 
147:                 POTLISTS(J1,2*J2-1,1) = ATOM1  ! Make an entry for each of these atoms in POTLISTS147:                 POTLISTS(J1,2*J2-1,1) = ATOM1  ! Make an entry for each of these atoms in POTLISTS
148:                 POTLISTS(J1,2*J2,1) = ATOM2148:                 POTLISTS(J1,2*J2,1) = ATOM2
149: 149: 
150:                 N_ATOM_PARTNERS(J1,2*J2-1)=1 ! ATOM1 partners with ATOM2...150:                 N_ATOM_PARTNERS(J1,2*J2-1)=1 ! ATOM1 partners with ATOM2...
151:                 POTLISTS(J1,2*J2-1,2) = ATOM2151:                 POTLISTS(J1,2*J2-1,2) = ATOM2
152:                 N_ATOM_PARTNERS(J1,2*J2)=0   ! ...but atom B doesn't partner with any (we've already included its interaction with A)152:                 N_ATOM_PARTNERS(J1,2*J2)=0   ! ...but atom B doesn't partner with any (we've already included its interaction with A)
213:                     POTLISTS(J1,J2,:N_ATOM_PARTNERS(J1,J2)+1)=ATOMLIST(J2:NATOM_BY_POT(J1))213:                     POTLISTS(J1,J2,:N_ATOM_PARTNERS(J1,J2)+1)=ATOMLIST(J2:NATOM_BY_POT(J1))
214:                 ENDDO214:                 ENDDO
215:             ENDIF215:             ENDIF
216: 216: 
217:             IF(DEBUG) THEN217:             IF(DEBUG) THEN
218:                 WRITE(*,*) "Potential:", POTTYPES(J1)218:                 WRITE(*,*) "Potential:", POTTYPES(J1)
219:                 WRITE(*,*) "N_ATOM_PARTNERS:"219:                 WRITE(*,*) "N_ATOM_PARTNERS:"
220:                 WRITE(*,*) N_ATOM_PARTNERS(J1,:NATOM_BY_POT(J1))220:                 WRITE(*,*) N_ATOM_PARTNERS(J1,:NATOM_BY_POT(J1))
221:                 WRITE(*,*) "POTLISTS:"221:                 WRITE(*,*) "POTLISTS:"
222:                 DO J2 = 1,NATOM_BY_POT(J1)222:                 DO J2 = 1,NATOM_BY_POT(J1)
223:                     IF (N_ATOM_PARTNERS(J1,J2).GT.0) THEN223:                     WRITE(*,*) "Atom number", J2, "in this potential"
224:                        WRITE(*,*) "Atom number", J2, "in this potential"224:                     WRITE(*,*) POTLISTS(J1,J2,N_ATOM_PARTNERS(J1,J2))
225:                        WRITE(*,*) POTLISTS(J1,J2,N_ATOM_PARTNERS(J1,J2)) 
226:                     ENDIF 
227:                 ENDDO225:                 ENDDO
228:             ENDIF226:             ENDIF
229: 227: 
230:         END SELECT228:         END SELECT
231:     ENDDO229:     ENDDO
232:     CLOSE(22)230:     CLOSE(22)
233: 231: 
234: END SUBROUTINE MULTIPOT_INITIALISE232: END SUBROUTINE MULTIPOT_INITIALISE
235: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!233: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
236: ! Calculate the potential energy and the gradient (if GTEST is set) and hessian (if STEST)234: ! Calculate the potential energy and the gradient (if GTEST is set) and hessian (if STEST)
268:                 CALL COMPUTE_PAIRWISE_POTENTIAL(X, G, ENERGY, GTEST, STEST, PAIRWISE_WCA, J1)266:                 CALL COMPUTE_PAIRWISE_POTENTIAL(X, G, ENERGY, GTEST, STEST, PAIRWISE_WCA, J1)
269:             CASE('HSPR')267:             CASE('HSPR')
270:                 ! Parameter is the equilibrium bond length, R0. Energy is returned in units of k_spr.268:                 ! Parameter is the equilibrium bond length, R0. Energy is returned in units of k_spr.
271:                 CALL COMPUTE_PAIRWISE_POTENTIAL(X, G, ENERGY, GTEST, STEST, HARMONIC_SPRINGS, J1)269:                 CALL COMPUTE_PAIRWISE_POTENTIAL(X, G, ENERGY, GTEST, STEST, HARMONIC_SPRINGS, J1)
272:             CASE DEFAULT270:             CASE DEFAULT
273:                 WRITE(*,*) "multipot> Error: unspecified potential"271:                 WRITE(*,*) "multipot> Error: unspecified potential"
274:                 STOP272:                 STOP
275:         END SELECT273:         END SELECT
276:     ENDDO274:     ENDDO
277: 275: 
278:  
279: END SUBROUTINE MULTIPOT_CALL276: END SUBROUTINE MULTIPOT_CALL
280: 277: 
281: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!278: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
282: !   Evaluate the energy(+gradient(+hessian)) for a set of atoms interacting according to this particular potential279: !   Evaluate the energy(+gradient(+hessian)) for a set of atoms interacting according to this particular potential
283: SUBROUTINE COMPUTE_ISOTROPIC_POTENTIAL(X, G, ENERGY, GTEST, STEST, POT, POTID, TMP_NATOMS)280: SUBROUTINE COMPUTE_ISOTROPIC_POTENTIAL(X, G, ENERGY, GTEST, STEST, POT, POTID, TMP_NATOMS)
284:     USE COMMONS, ONLY: NATOMS281:     USE COMMONS, ONLY: NATOMS
285:     USE MODHESS282:     USE MODHESS
286:     IMPLICIT NONE283:     IMPLICIT NONE
287:     DOUBLE PRECISION, INTENT(IN)    :: X(3*NATOMS)284:     DOUBLE PRECISION, INTENT(IN)    :: X(3*NATOMS)
288:     DOUBLE PRECISION, INTENT(INOUT) :: G(3*NATOMS)285:     DOUBLE PRECISION, INTENT(INOUT) :: G(3*NATOMS)


r30589/pairwise_potentials.f90 2016-06-13 18:30:06.972984611 +0100 r30588/pairwise_potentials.f90 2016-06-13 18:30:07.624993414 +0100
 23: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 23: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 24: ! Simple LJ potential, no cutoff. Atom radius sigma is given as a variable, but will often be set to 1. 24: ! Simple LJ potential, no cutoff. Atom radius sigma is given as a variable, but will often be set to 1.
 25: SUBROUTINE PAIRWISE_LJ(X1, X2, PARAMS, PG, PAIR_ENERGY, P_HESS, GTEST, STEST) 25: SUBROUTINE PAIRWISE_LJ(X1, X2, PARAMS, PG, PAIR_ENERGY, P_HESS, GTEST, STEST)
 26:     IMPLICIT NONE 26:     IMPLICIT NONE
 27:     DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3), PARAMS(10) 27:     DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3), PARAMS(10)
 28:     DOUBLE PRECISION, INTENT(OUT) :: PAIR_ENERGY 28:     DOUBLE PRECISION, INTENT(OUT) :: PAIR_ENERGY
 29:     DOUBLE PRECISION, INTENT(OUT) :: PG(6), P_HESS(6,6) 29:     DOUBLE PRECISION, INTENT(OUT) :: PG(6), P_HESS(6,6)
 30:     LOGICAL, INTENT(IN)           :: GTEST, STEST 30:     LOGICAL, INTENT(IN)           :: GTEST, STEST
 31:  31: 
 32:     DOUBLE PRECISION :: R2, R6, R8, R14, SIG, SIG6, SIG12  ! Various powers of the distance between the atoms, and the atom radius 32:     DOUBLE PRECISION :: R2, R6, R8, R14, SIG, SIG6, SIG12  ! Various powers of the distance between the atoms, and the atom radius
 33:     DOUBLE PRECISION :: G, F  ! G tensor and F tensor (see PAIRWISE_GRAD and PAIRWISE_HESSIAN, below) 33:         DOUBLE PRECISION :: G, F  ! G tensor and F tensor (see PAIRWISE_GRAD and PAIRWISE_HESSIAN, below)
 34:  34: 
 35:     SIG = PARAMS(1) 35:     SIG = PARAMS(1)
 36:     SIG6 = SIG**6 36:     SIG6 = SIG**6
 37:     SIG12 = SIG6**2 37:     SIG12 = SIG6**2
 38:  38: 
 39:     R2 = (X1(1)-X2(1))**2+(X1(2)-X2(2))**2+(X1(3)-X2(3))**2 39:     R2 = (X1(1)-X2(1))**2+(X1(2)-X2(2))**2+(X1(3)-X2(3))**2
 40:     R2 = 1.0D0/R2 40:     R2 = 1.0D0/R2
 41:     R6 = R2**3 41:     R6 = R2**3
 42:  42: 
 43:     PAIR_ENERGY=4.0D0*SIG6*R6*(SIG6*R6-1.0D0) 43:     PAIR_ENERGY=4.0D0*SIG6*R6*(SIG6*R6-1.0D0)
117: ! Hooke's Law spring potential117: ! Hooke's Law spring potential
118: ! The equilibrium bond length is given as PARAMS(1) (saved as R0). Energy is given in units of the spring constant.118: ! The equilibrium bond length is given as PARAMS(1) (saved as R0). Energy is given in units of the spring constant.
119: SUBROUTINE HARMONIC_SPRINGS(X1, X2, PARAMS, PG, PAIR_ENERGY, P_HESS, GTEST, STEST)119: SUBROUTINE HARMONIC_SPRINGS(X1, X2, PARAMS, PG, PAIR_ENERGY, P_HESS, GTEST, STEST)
120:     IMPLICIT NONE120:     IMPLICIT NONE
121:     DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3), PARAMS(10)121:     DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3), PARAMS(10)
122:     DOUBLE PRECISION, INTENT(OUT) :: PAIR_ENERGY122:     DOUBLE PRECISION, INTENT(OUT) :: PAIR_ENERGY
123:     DOUBLE PRECISION, INTENT(OUT) :: PG(6), P_HESS(6,6)123:     DOUBLE PRECISION, INTENT(OUT) :: PG(6), P_HESS(6,6)
124:     LOGICAL, INTENT(IN)           :: GTEST, STEST  ! GTEST is true when calculating the gradient as well as energy.124:     LOGICAL, INTENT(IN)           :: GTEST, STEST  ! GTEST is true when calculating the gradient as well as energy.
125:                                                        ! STEST is true when calculating the Hessian as well as the other two.125:                                                        ! STEST is true when calculating the Hessian as well as the other two.
126: 126: 
127:     DOUBLE PRECISION :: R0, R, R2, G, F127:         DOUBLE PRECISION :: R0, R, R2, G, F
128: 128: 
129:     R0 = PARAMS(1)129:     R0 = PARAMS(1)
130:     R2 = (X1(1)-X2(1))**2+(X1(2)-X2(2))**2+(X1(3)-X2(3))**2130:     R2 = (X1(1)-X2(1))**2+(X1(2)-X2(2))**2+(X1(3)-X2(3))**2
131:     R = SQRT(R2)131:     R = SQRT(R2)
132: 132: 
133:     PAIR_ENERGY = 0.5D0*(R-R0)**2133:     PAIR_ENERGY = 0.5D0*(R-R0)**2
134: 134: 
135:     IF(.NOT.GTEST) RETURN135:     IF(.NOT.GTEST) RETURN
136: 136: 
137:     F = R0/R137:     F = R0/R
138:     G = 1.0D0 - F138:     G = 1.0D0 - F
139:     R2 = 1.0D0/R2  ! PAIRWISE_HESSIAN requires the inverse square distance to be passed in. (cf. LJ/WCA potentials) 
140: 139: 
141:     CALL PAIRWISE_GRAD(X1,X2,G,PG)140:     CALL PAIRWISE_GRAD(X1,X2,G,PG)
142: 141: 
143:     IF (.NOT.STEST) RETURN142:     IF (.NOT.STEST) RETURN
144: 143: 
145:     CALL PAIRWISE_HESSIAN(X1,X2,G,F,R2,P_HESS)144:     CALL PAIRWISE_HESSIAN(X1,X2,G,F,R2,P_HESS)
146: 145: 
147:     RETURN146:     RETURN
148: 147: 
149: END SUBROUTINE HARMONIC_SPRINGS148: END SUBROUTINE HARMONIC_SPRINGS
171:     PG(4) = G*(X2(1)-X1(1))170:     PG(4) = G*(X2(1)-X1(1))
172:     PG(5) = G*(X2(2)-X1(2))171:     PG(5) = G*(X2(2)-X1(2))
173:     PG(6) = G*(X2(3)-X1(3))172:     PG(6) = G*(X2(3)-X1(3))
174: 173: 
175:     RETURN174:     RETURN
176: END SUBROUTINE PAIRWISE_GRAD175: END SUBROUTINE PAIRWISE_GRAD
177: 176: 
178: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!177: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
179: ! For all isotropic potentials the Hessian is calculated in the same way. We simply need to know the positions of the178: ! For all isotropic potentials the Hessian is calculated in the same way. We simply need to know the positions of the
180: ! two particles, the G tensor (see above) and the F tensor: F = r*d[(1/r)dU/dr]/dr179: ! two particles, the G tensor (see above) and the F tensor: F = r*d[(1/r)dU/dr]/dr
181: ! NOTE: R2 must be the INVERSE square distance! 
182: SUBROUTINE PAIRWISE_HESSIAN(X1,X2,G,F,R2,P_HESS)180: SUBROUTINE PAIRWISE_HESSIAN(X1,X2,G,F,R2,P_HESS)
183:     IMPLICIT NONE181:     IMPLICIT NONE
184:     DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3), G, F, R2182:     DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3), G, F, R2
185:     DOUBLE PRECISION, INTENT(OUT) :: P_HESS(6,6)183:     DOUBLE PRECISION, INTENT(OUT) :: P_HESS(6,6)
186:     INTEGER :: J1,J2184:     INTEGER :: J1,J2
187: 185: 
188: !       Compute the hessian. First are the entirely diagonal terms (6 of them)186: !       Compute the hessian. First are the entirely diagonal terms (6 of them)
189: !       These correspond to 2nd derivatives wrt the same coordinate twice187: !       These correspond to 2nd derivatives wrt the same coordinate twice
190:     DO J1=1,3188:     DO J1=1,3
191:         P_HESS(J1,J1) = F*R2*(X1(J1)-X2(J1))**2 + G189:         P_HESS(J1,J1) = F*R2*(X1(J1)-X2(J1))**2 + G


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0