hdiff output

r31189/isotropic_potentials.f90 2016-09-23 11:30:21.834874783 +0100 r31188/isotropic_potentials.f90 2016-09-23 11:30:22.602885020 +0100
253:     CALL ISOTROPIC_GRAD(N, X, G, TMP_G)253:     CALL ISOTROPIC_GRAD(N, X, G, TMP_G)
254: 254: 
255:     IF (.NOT.STEST) RETURN  ! It is assumed we will never need the Hessian without also needing the gradient.255:     IF (.NOT.STEST) RETURN  ! It is assumed we will never need the Hessian without also needing the gradient.
256: 256: 
257:     CALL ISOTROPIC_HESSIAN(N, X, G, F, R2, TMP_HESS)257:     CALL ISOTROPIC_HESSIAN(N, X, G, F, R2, TMP_HESS)
258: 258: 
259:     RETURN259:     RETURN
260: 260: 
261: END SUBROUTINE ISO_WCA261: END SUBROUTINE ISO_WCA
262: 262: 
263: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
264: ! Coulomb potential (no cutoff - use for nonperiodic potentials only) 
265: ! This subroutine is called directly from multipot.f90, so that we can avoid calculating unnecessary interactions 
266: ! (see below). The same approach is used for the "EXCLUDE" potentials. 
267: SUBROUTINE ISO_COULOMB(X, POTLIST, N_ATOM_PARTNERS, Q, POTSCALE, PARAMS, ENERGY, GRAD, GTEST, STEST) 
268:  
269:     ! Because we're passing in the entire coordinates array here, we need to know the total number of atoms 
270:     USE COMMONS, ONLY: NATOMS 
271:  
272:     IMPLICIT NONE 
273:  
274:     DOUBLE PRECISION, INTENT(IN)  :: X(3*NATOMS)          ! The full atomistic coordinate array 
275:     INTEGER, INTENT(IN)           :: POTLIST(:,:)         ! An array containing the interaction partners for each atom 
276:     INTEGER, INTENT(IN)           :: N_ATOM_PARTNERS(:)   ! An array containing the number of interaction partners for each atom 
277:     DOUBLE PRECISION, INTENT(IN)  :: Q (:)                ! An array containing the atom charges 
278:     DOUBLE PRECISION, INTENT(IN)  :: POTSCALE             ! The energy unit for this potential, which multiplies E, G and Hess. 
279:     DOUBLE PRECISION, INTENT(IN)  :: PARAMS(10)           ! Maximum number of parameters is hardcoded here 
280:     DOUBLE PRECISION, INTENT(OUT) :: ENERGY               ! The energy of the configuration 
281:     DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOMS)       ! The energy gradient 
282:     LOGICAL, INTENT(IN)           :: GTEST, STEST         ! Flags to specify whether the gradient should be calculated (GTEST) and 
283:                                                           ! whether the Hessian should be calculated (STEST) 
284:  
285:     ! Various powers of the distance between the atoms, and the atom radius 
286:     DOUBLE PRECISION :: R, R2(NATOMS,NATOMS), R3 
287:     DOUBLE PRECISION :: G(NATOMS,NATOMS), F(NATOMS,NATOMS)  ! G tensor and F tensor (see ISOTROPIC_GRAD and ISOTROPIC_HESSIAN, below) 
288:     DOUBLE PRECISION :: TMP_ENERGY 
289:     INTEGER :: J1, J2, J3, J4, J5, J6 
290:  
291:     TMP_ENERGY=0.0D0 
292:                  
293:     ! The arrangement of these IF statements seems slightly odd, but it's been done to make sure we only need to set up one pair 
294:     ! of DO loops for each call. 
295:     IF (GTEST) THEN 
296:         IF (STEST) THEN  ! Gradient + Hessian 
297:  
298:             ! Rather than looping over all pairs of atoms, we loop over all interacting pairs as specified by POTLIST. 
299:             ! The first entry in each element of POTLIST contains the index of a principal atom. All other entries contain 
300:             ! the interaction partners for this principal atom. 
301:             DO J6=1,UBOUND(POTLIST,1) 
302:                 J1=POTLIST(J6,1)  ! Get the actual atomic index for this atom 
303:                 J3=3*J1 
304:                 R2(J1,J1)=0.0D0 
305:                 G(J1,J1)=0.0D0 
306:                 F(J1,J1)=0.0D0 
307:  
308:                 ! Loop over the partners for atom J1 
309:                 DO J5=2,N_ATOM_PARTNERS(J6)+1  ! Have to add the +1 because the first element in POTLIST(J1,:) 
310:                                                ! is the atom index J1, instead of a neighbour index. 
311:  
312:                     J2=POTLIST(J6,J5)  ! Get the actual atomic index for this atom 
313:                     ! For historical reasons, the order of the variables J2, J4, J5 is odd. 
314:                     J4=3*J2 
315:                     
316:                     R2(J2,J1)=(X(J3-2)-X(J4-2))**2+(X(J3-1)-X(J4-1))**2+(X(J3)-X(J4))**2 
317:                     R2(J2,J1)=1.0D0/R2(J2,J1)  ! From now on, both R2 and R are inverse lengths 
318:                     R=SQRT(R2(J2,J1)) 
319:  
320:                     TMP_ENERGY=TMP_ENERGY+Q(J1)*Q(J2)*R 
321:  
322: !                    write(*,*) "E+G+H. Pair, Dist, Energy:", J1, J2, 1.0D0/R2(J2,J1), SIG6*R6*(SIG6*R6-1.0D0) + 0.25D0 
323:                     ! Set up storage arrays to use in the gradient- and hessian-calculating routines 
324:                     R2(J1,J2) = R2(J2,J1) 
325:                     R3 = R*R2(J2,J1) 
326:  
327:                     G(J2,J1)=-Q(J1)*Q(J2)*R3 
328:                     G(J1,J2)=G(J2,J1) 
329:  
330:                     F(J2,J1)=-3*G(J2,J1) 
331:                     F(J1,J2)=F(J2,J1) 
332:  
333:                 ENDDO 
334:             ENDDO 
335:         ELSE             ! Energy + Gradient only 
336:             DO J6=1,UBOUND(POTLIST,1) 
337:                 J1=POTLIST(J6,1) 
338:                 J3=3*J1 
339:                 G(J1,J1)=0.0D0 
340:  
341:                 DO J5=2,N_ATOM_PARTNERS(J6)+1 
342:                     
343:                     J2=POTLIST(J6,J5)  ! For historical reasons, the order of the variables J2, J4, J5 is odd. 
344:                     J4=3*J2 
345:  
346:                     R=SQRT((X(J3-2)-X(J4-2))**2+(X(J3-1)-X(J4-1))**2+(X(J3)-X(J4))**2) 
347:                     R=1.0D0/R 
348:                     R3=R**3 
349:  
350:                     TMP_ENERGY=TMP_ENERGY+Q(J1)*Q(J2)*R 
351: !                    write(*,*) "E+G. Pair, Dist, Energy:", J1, J2, (X(J3-2)-X(J4-2))**2+(X(J3-1)-X(J4-1))**2+(X(J3)-X(J4))**2, SIG6*R6*(SIG6*R6-1.0D0) + 0.25D0 
352:                     G(J2,J1)=-Q(J1)*Q(J2)*R3 
353:                     G(J1,J2)=G(J2,J1) 
354:                 ENDDO 
355:             ENDDO 
356:         ENDIF 
357:     ELSE                ! Energy only 
358:         DO J6=1,UBOUND(POTLIST,1) 
359:             J1=POTLIST(J6,1) 
360:             J3=3*J1 
361:  
362:             DO J5=2,N_ATOM_PARTNERS(J6)+1 
363:  
364:                 J2=POTLIST(J6,J5)  ! For historical reasons, the order of the variables J2, J4, J5 is odd. 
365:                 J4=3*J2 
366:  
367:                 R=SQRT((X(J3-2)-X(J4-2))**2+(X(J3-1)-X(J4-1))**2+(X(J3)-X(J4))**2) 
368:                 R=1.0D0/R 
369:  
370:                 TMP_ENERGY=TMP_ENERGY+Q(J1)*Q(J2)*R 
371: !                write(*,*) "E. Pair, Dist, Energy:", J1, J2, R2(J2,J1), SIG6*R6*(SIG6*R6-1.0D0) + 0.25D0 
372:             ENDDO 
373:         ENDDO 
374:  
375:     ENDIF 
376:  
377:     ENERGY = ENERGY + POTSCALE*TMP_ENERGY 
378:  
379:     IF (.NOT.GTEST) RETURN 
380:  
381:     ! G should already be set 
382:     CALL EXCLUDE_ISOTROPIC_GRAD(POTLIST, N_ATOM_PARTNERS, X, G, GRAD, POTSCALE) 
383:  
384:     IF (.NOT.STEST) RETURN  ! It is assumed we will never need the Hessian without also needing the gradient. 
385:  
386:     ! This is currently a very inefficient implementation! I haven't implemented the use of POTLISTS yet. 
387:     CALL EXCLUDE_ISOTROPIC_HESSIAN(POTLIST, N_ATOM_PARTNERS, X, G, F, R2, POTSCALE) 
388:  
389:     RETURN 
390:  
391: END SUBROUTINE ISO_COULOMB 
392:  
393:  
394:  
395:  
396: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!263: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
397: ! LJ potential with support for excluded interaction sets. 264: ! LJ potential with support for excluded interaction sets. 
398: ! Atom radius sigma is given as a variable, but will often be set to 1.265: ! Atom radius sigma is given as a variable, but will often be set to 1.
399: SUBROUTINE EXCLUDE_ISO_LJ(X, POTLIST, N_ATOM_PARTNERS, POTSCALE, PARAMS, ENERGY, GRAD, GTEST, STEST)266: SUBROUTINE EXCLUDE_ISO_LJ(X, POTLIST, N_ATOM_PARTNERS, POTSCALE, PARAMS, ENERGY, GRAD, GTEST, STEST)
400: 267: 
401:     ! Because we're passing in the entire coordinates array here, we need to know the total number of atoms268:     ! Because we're passing in the entire coordinates array here, we need to know the total number of atoms
402:     USE COMMONS, ONLY: NATOMS269:     USE COMMONS, ONLY: NATOMS
403: !    USE MODHESS270: !    USE MODHESS
404: 271: 
405:     IMPLICIT NONE272:     IMPLICIT NONE
731: 598: 
732:     RETURN599:     RETURN
733: END SUBROUTINE ISOTROPIC_GRAD600: END SUBROUTINE ISOTROPIC_GRAD
734: 601: 
735: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!602: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
736: ! For all isotropic potentials the Hessian is calculated in the same way. We simply need to know the positions of the603: ! For all isotropic potentials the Hessian is calculated in the same way. We simply need to know the positions of the
737: ! two particles, the G tensor (see above) and the F tensor: F = r*d[(1/r)dU/dr]/dr604: ! two particles, the G tensor (see above) and the F tensor: F = r*d[(1/r)dU/dr]/dr
738: SUBROUTINE ISOTROPIC_HESSIAN(N, X, G, F, R2, TMP_HESS)605: SUBROUTINE ISOTROPIC_HESSIAN(N, X, G, F, R2, TMP_HESS)
739:     IMPLICIT NONE606:     IMPLICIT NONE
740:     INTEGER, INTENT(IN)           :: N  ! Number of atoms607:     INTEGER, INTENT(IN)           :: N  ! Number of atoms
741:     DOUBLE PRECISION, INTENT(IN)  :: X(3*N), G(N,N), F(N,N), R2(N,N) ! NB: R2 must be an inverse square distance!!!!!608:     DOUBLE PRECISION, INTENT(IN)  :: X(3*N), G(N,N), F(N,N), R2(N,N)
742:                                                                      ! This is a consequence of the way I've defined F. 
743:     DOUBLE PRECISION, INTENT(OUT) :: TMP_HESS(3*N,3*N)609:     DOUBLE PRECISION, INTENT(OUT) :: TMP_HESS(3*N,3*N)
744:     INTEGER :: J1, J2, J3, J4, J5, J6610:     INTEGER :: J1, J2, J3, J4, J5, J6
745:     DOUBLE PRECISION :: DUMMY611:     DOUBLE PRECISION :: DUMMY
746: 612: 
747: !       Compute the hessian. First are the entirely diagonal terms (3N of them)613: !       Compute the hessian. First are the entirely diagonal terms (3N of them)
748: !       These correspond to 2nd derivatives wrt the same coordinate twice614: !       These correspond to 2nd derivatives wrt the same coordinate twice
749:     DO J1=1,N615:     DO J1=1,N
750:         DO J2=1,3616:         DO J2=1,3
751:             J3=3*(J1-1)+J2617:             J3=3*(J1-1)+J2
752:             DUMMY=0.0D0618:             DUMMY=0.0D0
772:             ENDDO638:             ENDDO
773:         ENDDO639:         ENDDO
774:     ENDDO640:     ENDDO
775: 641: 
776: !  Case III, different atoms, same cartesian coordinate.642: !  Case III, different atoms, same cartesian coordinate.
777: 643: 
778:     DO J1=1,N644:     DO J1=1,N
779:         DO J2=1,3645:         DO J2=1,3
780:             J3=3*(J1-1)+J2646:             J3=3*(J1-1)+J2
781:             DO J4=J1+1,N647:             DO J4=J1+1,N
782:                ! Is it ok that I'm not initialising to 0 and using the TMP_HESS construction here? 
783:                ! I need to check this! 
784:                 TMP_HESS(3*(J4-1)+J2,J3)=-F(J4,J1)*R2(J4,J1)*(X(J3)-X(3*(J4-1)+J2))**2-G(J4,J1)648:                 TMP_HESS(3*(J4-1)+J2,J3)=-F(J4,J1)*R2(J4,J1)*(X(J3)-X(3*(J4-1)+J2))**2-G(J4,J1)
785:             ENDDO649:             ENDDO
786:         ENDDO650:         ENDDO
787:     ENDDO651:     ENDDO
788: 652: 
789: !  Case IV: different atoms and different cartesian coordinates.653: !  Case IV: different atoms and different cartesian coordinates.
790: 654: 
791:     DO J1=1,N655:     DO J1=1,N
792:         DO J2=1,3656:         DO J2=1,3
793:             J3=3*(J1-1)+J2657:             J3=3*(J1-1)+J2
809:         ENDDO673:         ENDDO
810:     ENDDO674:     ENDDO
811: 675: 
812:     RETURN676:     RETURN
813: END SUBROUTINE ISOTROPIC_HESSIAN677: END SUBROUTINE ISOTROPIC_HESSIAN
814: 678: 
815: 679: 
816: ! For all isotropic potentials the gradient is calculated in the same way. We simply need to know the positions of the680: ! For all isotropic potentials the gradient is calculated in the same way. We simply need to know the positions of the
817: ! two particles, and the G tensor: G = (1/r)dU/dr for an isotropic potential U(r)681: ! two particles, and the G tensor: G = (1/r)dU/dr for an isotropic potential U(r)
818: SUBROUTINE EXCLUDE_ISOTROPIC_GRAD(POTLIST, N_ATOM_PARTNERS, X, G, GRAD, POTSCALE)682: SUBROUTINE EXCLUDE_ISOTROPIC_GRAD(POTLIST, N_ATOM_PARTNERS, X, G, GRAD, POTSCALE)
 683: 
819:     IMPLICIT NONE684:     IMPLICIT NONE
820:     INTEGER, INTENT(IN)           :: POTLIST(:,:)685:     INTEGER, INTENT(IN)           :: POTLIST(:,:)
821:     INTEGER, INTENT(IN)           :: N_ATOM_PARTNERS(:)  ! Number of atoms686:     INTEGER, INTENT(IN)           :: N_ATOM_PARTNERS(:)  ! Number of atoms
822:     DOUBLE PRECISION, INTENT(IN)  :: X(:), G(:,:)687:     DOUBLE PRECISION, INTENT(IN)  :: X(:), G(:,:)
823:     DOUBLE PRECISION, INTENT(OUT) :: GRAD(:)  ! This is the real full gradient, so we don't want to zero it!688:     DOUBLE PRECISION, INTENT(OUT) :: GRAD(:)  ! This is the real full gradient, so we don't want to zero it!
824:     DOUBLE PRECISION, INTENT(IN)  :: POTSCALE689:     DOUBLE PRECISION, INTENT(IN)  :: POTSCALE
825:     INTEGER :: J1, J2, J3, J4, J5, J6690:     INTEGER :: J1, J2, J3, J4, J5, J6
826:     DOUBLE PRECISION :: DUMMYX, DUMMYY, DUMMYZ, XMUL691:     DOUBLE PRECISION :: DUMMYX, DUMMYY, DUMMYZ, XMUL
827: 692: 
828:     DO J6=1,UBOUND(POTLIST,1)  ! The atom for which we are calculating the gradient693:     DO J6=1,UBOUND(POTLIST,1)  ! The atom for which we are calculating the gradient


r31189/multipot.f90 2016-09-23 11:30:22.090878196 +0100 r31188/multipot.f90 2016-09-23 11:30:22.854888378 +0100
 33:     INTEGER, ALLOCATABLE :: POTLISTS(:,:,:)         ! Holds the indices of all atoms belonging to each potential type, and 33:     INTEGER, ALLOCATABLE :: POTLISTS(:,:,:)         ! Holds the indices of all atoms belonging to each potential type, and
 34:                                                     ! the indices of the partners which interact with each atom 34:                                                     ! the indices of the partners which interact with each atom
 35:                                                     ! POTLISTS(m,n,1) contains the index of the nth atom which has an attraction of 35:                                                     ! POTLISTS(m,n,1) contains the index of the nth atom which has an attraction of
 36:                                                     ! type POTTYPES(m), and POTLISTS(m,n,2:N_ATOM_PARTNERS(m,n)) are the indices of 36:                                                     ! type POTTYPES(m), and POTLISTS(m,n,2:N_ATOM_PARTNERS(m,n)) are the indices of
 37:                                                     ! its partners. 37:                                                     ! its partners.
 38:  38: 
 39:     INTEGER, ALLOCATABLE :: DEGS_BY_POT(:,:)        ! For use by isotropic potentials, where all atoms interact with all other atoms. 39:     INTEGER, ALLOCATABLE :: DEGS_BY_POT(:,:)        ! For use by isotropic potentials, where all atoms interact with all other atoms.
 40:                                                     ! Contains a list of the degrees of freedom belonging to all the atoms using this 40:                                                     ! Contains a list of the degrees of freedom belonging to all the atoms using this
 41:                                                     ! potential. 41:                                                     ! potential.
 42:  42: 
 43:     DOUBLE PRECISION, ALLOCATABLE :: CHRG(:)        ! Charges for individual atoms. Used by ISO_COULOMB. 
 44:  
 45:  43: 
 46: CONTAINS 44: CONTAINS
 47: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 45: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 48:     ! Parse the input file which is required to specify the different types of interaction we have in the system 46:     ! Parse the input file which is required to specify the different types of interaction we have in the system
 49: SUBROUTINE MULTIPOT_INITIALISE 47: SUBROUTINE MULTIPOT_INITIALISE
 50:  48: 
 51:     USE COMMONS, ONLY: NATOMS 49:     USE COMMONS, ONLY: NATOMS
 52:     USE KEY, ONLY: BULKT 
 53:     USE GENRIGID, ONLY: RIGIDINIT, RB_BY_ATOM 50:     USE GENRIGID, ONLY: RIGIDINIT, RB_BY_ATOM
 54:  51: 
 55:     IMPLICIT NONE 52:     IMPLICIT NONE
 56:  53: 
 57:     INTEGER NPARAMS, ATOM1, ATOM2 54:     INTEGER NPARAMS, ATOM1, ATOM2
 58:     INTEGER :: ATOMLIST(NATOMS) 55:     INTEGER :: ATOMLIST(NATOMS)
 59:     CHARACTER(LEN=1000) :: DUMMYCHAR 56:     CHARACTER(LEN=1000) :: DUMMYCHAR
 60:     LOGICAL :: END 57:     LOGICAL :: END
 61:     INTEGER :: J1, J2, J3, J4, iostatus, COUNTER 58:     INTEGER :: J1, J2, J3, J4, iostatus, COUNTER
 62:  59: 
168:                 WRITE(*,*) "Potential:", POTTYPES(J1)165:                 WRITE(*,*) "Potential:", POTTYPES(J1)
169:                 WRITE(*,*) "N_ATOM_PARTNERS:"166:                 WRITE(*,*) "N_ATOM_PARTNERS:"
170:                 WRITE(*,*) N_ATOM_PARTNERS(J1,:NATOM_BY_POT(J1))167:                 WRITE(*,*) N_ATOM_PARTNERS(J1,:NATOM_BY_POT(J1))
171:                 WRITE(*,*) "POTLISTS:"168:                 WRITE(*,*) "POTLISTS:"
172:                 DO J2 = 1,NATOM_BY_POT(J1)169:                 DO J2 = 1,NATOM_BY_POT(J1)
173:                     WRITE(*,*) "Atom number", J2, "in this potential"170:                     WRITE(*,*) "Atom number", J2, "in this potential"
174:                     WRITE(*,*) POTLISTS(J1,J2,:N_ATOM_PARTNERS(J1,J2))171:                     WRITE(*,*) POTLISTS(J1,J2,:N_ATOM_PARTNERS(J1,J2))
175:                 ENDDO172:                 ENDDO
176:             ENDIF173:             ENDIF
177: 174: 
178:         CASE('ILJ','IWCA') 
179:         ! For "iso_" potentials. Every specified atom with this potential type interacts with all the others.175:         ! For "iso_" potentials. Every specified atom with this potential type interacts with all the others.
180:         ! The input format is simple: list all the atoms using this potential on a single line, separated by spaces.176:         ! The input format is simple: list all the atoms using this potential on a single line, separated by spaces.
181:         ! They don't have to be in index order, but everything will make more sense if they are!177:         ! They don't have to be in index order, but everything will make more sense if they are!
182:         ! There is currently no facility to exclude the interactions within a rigid body. Use EWCA/ELJ if you need that.178:         ! There is currently no facility to exclude the interactions within a rigid body. Use EWCA/ELJ if you need that.
183: 179:         CASE('ILJ','IWCA')
184:             READ(22,*) ATOMLIST(:NATOM_BY_POT(J1))180:             READ(22,*) ATOMLIST(:NATOM_BY_POT(J1))
185: 181: 
186:             ! All we need to save for this type of potential is the list of whole-system degrees of freedom on which182:             ! All we need to save for this type of potential is the list of whole-system degrees of freedom on which
187:             ! the potential will depend.183:             ! the potential will depend.
188:             DO J2=1,NATOM_BY_POT(J1)184:             DO J2=1,NATOM_BY_POT(J1)
189:                 J3 = ATOMLIST(J2)185:                 J3 = ATOMLIST(J2)
190:                 DEGS_BY_POT(J1,3*J2-2) = 3*J3-2186:                 DEGS_BY_POT(J1,3*J2-2) = 3*J3-2
191:                 DEGS_BY_POT(J1,3*J2-1) = 3*J3-1187:                 DEGS_BY_POT(J1,3*J2-1) = 3*J3-1
192:                 DEGS_BY_POT(J1,3*J2)   = 3*J3188:                 DEGS_BY_POT(J1,3*J2)   = 3*J3
193:             ENDDO189:             ENDDO
194: 190: 
195:             IF(DEBUG) THEN191:             IF(DEBUG) THEN
196:                 write(*,*) "Potential:", POTTYPES(J1)192:                 write(*,*) "Potential:", POTTYPES(J1)
197:                 write(*,*) "DEGS_BY_POT:"193:                 write(*,*) "DEGS_BY_POT:"
198:                 write(*,*) DEGS_BY_POT(J1,:NATOM_BY_POT(J1)*3)194:                 write(*,*) DEGS_BY_POT(J1,:NATOM_BY_POT(J1)*3)
199:             ENDIF195:             ENDIF
200: 196: 
201:         CASE('ICOU') 
202:         ! Coulomb potential with no truncation or smoothing. 
203:         ! Input format: list each charged atom on its own line, followed by its charge. 
204:         ! If RIGIDINIT is set, interactions within a rigid body are not calculated. 
205:             
206:            IF(BULKT) THEN 
207:               WRITE(*,*) "multipot> ERROR: Isotropic Coulomb potential has been specified with periodic boundaries." 
208:               WRITE(*,*) "PBCs are not implemented for this potential." 
209:               STOP 
210:            ENDIF 
211:  
212:            ! Initialise charge array - start by assuming all atoms have 0 charge. 
213:            ALLOCATE(CHRG(NATOMS)) 
214:            CHRG(:) = 0.0D0 
215:            N_ATOM_PARTNERS(J1,:)=0 
216:  
217:            DO J2=1,NATOM_BY_POT(J1) 
218:               READ(22,*) POTLISTS(J1,J2,1), CHRG(POTLISTS(J1,J2,1)) 
219:            ENDDO 
220:  
221:            IF (RIGIDINIT) THEN  ! Need to avoid adding interactions within a rigid body. 
222:  
223:               IF (.NOT. ALLOCATED(RB_BY_ATOM)) THEN 
224:                  WRITE(*,*) "multipot> Error: RB_BY_ATOM is not allocated by the time MULTIPOT_INITIALISE is called." 
225:                  WRITE(*,*) "multipot> The RIGIDINIT keyword must precede the MULTIPOT keyword in your odata file." 
226:                  STOP 98 
227:               ENDIF 
228:  
229:               DO J2=1,NATOM_BY_POT(J1) 
230:                  ATOM1 = POTLISTS(J1,J2,1) 
231:                  DO J3=J2+1,NATOM_BY_POT(J1) ! The lower limit here ensures we don't add a pair twice, and no atom  
232:                                              ! interacts with itself 
233:                     ATOM2 = POTLISTS(J1,J3,1) 
234:                     IF(RB_BY_ATOM(ATOM1) .NE. RB_BY_ATOM(ATOM2)) THEN ! The atoms are in different rigid bodies 
235:                        N_ATOM_PARTNERS(J1,J2) = N_ATOM_PARTNERS(J1,J2) + 1 
236:                        POTLISTS(J1,J2,1+N_ATOM_PARTNERS(J1,J2)) = ATOM2 
237:                     ENDIF 
238:                  ENDDO 
239:               ENDDO 
240:  
241:            ELSE ! Then every atom partners with every other atom in this potential exactly once. 
242:               DO J2=1,NATOM_BY_POT(J1) 
243:                  DO J3=J2+1,NATOM_BY_POT(J1) ! The lower limit here ensures we don't add a pair twice. 
244:                     N_ATOM_PARTNERS(J1,J2) = N_ATOM_PARTNERS(J1,J2) + 1 
245:                     POTLISTS(J1,J2,1+N_ATOM_PARTNERS(J1,J2)) = POTLISTS(J1,J3,1) 
246:                  ENDDO 
247:               ENDDO 
248:            ENDIF 
249:  
250:            IF(DEBUG) THEN 
251:               WRITE(*,*) "Potential:", POTTYPES(J1) 
252:               WRITE(*,*) "Number of atoms:", NATOM_BY_POT(J1) 
253:               WRITE(*,*) "N_ATOM_PARTNERS:" 
254:               WRITE(*,*) N_ATOM_PARTNERS(J1,:NATOM_BY_POT(J1)) 
255:               WRITE(*,*) "POTLISTS:" 
256:               DO J2 = 1,NATOM_BY_POT(J1) 
257:                  WRITE(*,*) "Atom number", J2, "in this potential" 
258:                  WRITE(*,*) POTLISTS(J1,J2,:1+N_ATOM_PARTNERS(J1,J2)) 
259:               ENDDO 
260:            ENDIF 
261:  
262:         CASE('EWCA', 'ELJ')197:         CASE('EWCA', 'ELJ')
263: 198: 
264:         ! An isotropic potential with the option to switch off particular interactions.199:         ! An isotropic potential with the option to switch off particular interactions.
265:         !200:         !
266:         ! Specify the allowed interactions as normal for isotropic potentials: with a list of atom indices on a single line.201:         ! Specify the allowed interactions as normal for isotropic potentials: with a list of atom indices on a single line.
267:         ! If RIGIDINIT is set, then interactions within a rigid body are not calculated.202:         ! If RIGIDINIT is set, then interactions within a rigid body are not calculated.
268:         ! Additional excluded interactions are specified by appending the following to this multipotconfig entry:203:         ! Additional excluded interactions are specified by appending the following to this multipotconfig entry:
269:         !204:         !
270:         ! EXCLUDE nlines max_length205:         ! EXCLUDE nlines max_length
271:         ! line_len_1206:         ! line_len_1
426:                 ! Parameter is the equilibrium bond length, R0. Energy is returned in units of k_spr.361:                 ! Parameter is the equilibrium bond length, R0. Energy is returned in units of k_spr.
427:                 CALL COMPUTE_PAIRWISE_POTENTIAL(X, G, ENERGY, GTEST, STEST, HARMONIC_SPRINGS, J1)362:                 CALL COMPUTE_PAIRWISE_POTENTIAL(X, G, ENERGY, GTEST, STEST, HARMONIC_SPRINGS, J1)
428:             ! For exclusion potentials, we must also pass in a list specifying the number of interacting partners possessed363:             ! For exclusion potentials, we must also pass in a list specifying the number of interacting partners possessed
429:             ! by each atom, and a list specifying which atoms make up these partners (from POTLISTS)364:             ! by each atom, and a list specifying which atoms make up these partners (from POTLISTS)
430:             CASE('ELJ')365:             CASE('ELJ')
431:                 ! Only one parameter for this potential: the particle radius, sigma366:                 ! Only one parameter for this potential: the particle radius, sigma
432:                 CALL EXCLUDE_ISO_LJ(X, POTLISTS(J1,:THISN,:), N_ATOM_PARTNERS(J1,:THISN), POTSCALES(J1), THESEPARAMS, ENERGY, G, GTEST, STEST)367:                 CALL EXCLUDE_ISO_LJ(X, POTLISTS(J1,:THISN,:), N_ATOM_PARTNERS(J1,:THISN), POTSCALES(J1), THESEPARAMS, ENERGY, G, GTEST, STEST)
433:             CASE('EWCA')368:             CASE('EWCA')
434:                  ! Only one parameter for this potential: the particle radius, sigma369:                  ! Only one parameter for this potential: the particle radius, sigma
435:                 CALL EXCLUDE_ISO_WCA(X, POTLISTS(J1,:THISN,:), N_ATOM_PARTNERS(J1,:THISN), POTSCALES(J1), THESEPARAMS, ENERGY, G, GTEST, STEST)370:                 CALL EXCLUDE_ISO_WCA(X, POTLISTS(J1,:THISN,:), N_ATOM_PARTNERS(J1,:THISN), POTSCALES(J1), THESEPARAMS, ENERGY, G, GTEST, STEST)
436:             CASE('ICOU') 
437:                 CALL ISO_COULOMB(X, POTLISTS(J1,:THISN,:), N_ATOM_PARTNERS(J1,:THISN), CHRG, POTSCALES(J1), THESEPARAMS, ENERGY, G, GTEST, STEST) 
438:             CASE DEFAULT371:             CASE DEFAULT
439:                 ! We shouldn't every get here, unless you have added a new type of potential to MULTIPOT_INITIALISE and forgotten372:                 ! We shouldn't every get here, unless you have added a new type of potential to MULTIPOT_INITIALISE and forgotten
440:                 ! to add it to this SELECT CASE.373:                 ! to add it to this SELECT CASE.
441:                 WRITE(*,*) "multipot> Error: unspecified potential"374:                 WRITE(*,*) "multipot> Error: unspecified potential"
442:                 STOP375:                 STOP
443:         END SELECT376:         END SELECT
444:     ENDDO377:     ENDDO
445: 378: 
446: END SUBROUTINE MULTIPOT_CALL379: END SUBROUTINE MULTIPOT_CALL
447: 380: 


r31189/pairwise_potentials.f90 2016-09-23 11:30:22.338881501 +0100 r31188/pairwise_potentials.f90 2016-09-23 11:30:23.102891685 +0100
115: 115: 
116: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!116: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0