hdiff output

r30623/isotropic_potentials.f90 2016-07-06 15:37:53.757093078 +0100 r30622/isotropic_potentials.f90 2016-07-06 15:37:55.941122641 +0100
  7: !            INTEGER, INTENT(IN)           :: TMP_NATOMS  7: !            INTEGER, INTENT(IN)           :: TMP_NATOMS
  8: !            DOUBLE PRECISION, INTENT(IN)  :: X(3*TMP_NATOMS)  8: !            DOUBLE PRECISION, INTENT(IN)  :: X(3*TMP_NATOMS)
  9: !            DOUBLE PRECISION, INTENT(IN)  :: PARAMS(10)      ! Maximum number of parameters is hardcoded here  9: !            DOUBLE PRECISION, INTENT(IN)  :: PARAMS(10)      ! Maximum number of parameters is hardcoded here
 10: !                                                             ! Each potential will use a subset of the elements of PARAMS 10: !                                                             ! Each potential will use a subset of the elements of PARAMS
 11: !            DOUBLE PRECISION, INTENT(OUT) :: TMP_ENERGY 11: !            DOUBLE PRECISION, INTENT(OUT) :: TMP_ENERGY
 12: !            DOUBLE PRECISION, INTENT(OUT) :: TMP_G(3*TMP_NATOMS), TMP_HESS(3*TMP_NATOMS,3*TMP_NATOMS) 12: !            DOUBLE PRECISION, INTENT(OUT) :: TMP_G(3*TMP_NATOMS), TMP_HESS(3*TMP_NATOMS,3*TMP_NATOMS)
 13: !            LOGICAL, INTENT(IN)           :: GTEST, STEST    ! GTEST is true when calculating the gradient as well as energy. 13: !            LOGICAL, INTENT(IN)           :: GTEST, STEST    ! GTEST is true when calculating the gradient as well as energy.
 14: !                                                             ! STEST is true when calculating the Hessian as well as the other two. 14: !                                                             ! STEST is true when calculating the Hessian as well as the other two.
 15: !        END SUBROUTINE POT 15: !        END SUBROUTINE POT
 16:  16: 
 17: ! An exception to this rule is the EXCLUDE potentials, which need additional information passed in. They are handled differently 
 18: ! (see multipot.f90) 
 19:  
 20: MODULE ISOTROPIC_POTENTIALS 17: MODULE ISOTROPIC_POTENTIALS
 21:  18: 
 22: IMPLICIT NONE 19: IMPLICIT NONE
 23:  20: 
 24: DOUBLE PRECISION, PARAMETER :: WCA_CUT=1.2599210498948731647672106072782283505702514647015079 ! = 2^(1/3) 21: DOUBLE PRECISION, PARAMETER :: WCA_CUT=1.2599210498948731647672106072782283505702514647015079 ! = 2^(1/3)
 25:  22: 
 26: CONTAINS 23: CONTAINS
 27:  24: 
 28: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 25: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 29: ! Simple LJ potential, no cutoff. Atom radius sigma is given as a variable, but will often be set to 1. 26: ! Simple LJ potential, no cutoff. Atom radius sigma is given as a variable, but will often be set to 1.
254: 251: 
255:     IF (.NOT.STEST) RETURN  ! It is assumed we will never need the Hessian without also needing the gradient.252:     IF (.NOT.STEST) RETURN  ! It is assumed we will never need the Hessian without also needing the gradient.
256: 253: 
257:     CALL ISOTROPIC_HESSIAN(N, X, G, F, R2, TMP_HESS)254:     CALL ISOTROPIC_HESSIAN(N, X, G, F, R2, TMP_HESS)
258: 255: 
259:     RETURN256:     RETURN
260: 257: 
261: END SUBROUTINE ISO_WCA258: END SUBROUTINE ISO_WCA
262: 259: 
263: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!260: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
264: ! LJ potential with support for excluded interaction sets.  
265: ! Atom radius sigma is given as a variable, but will often be set to 1. 
266: SUBROUTINE EXCLUDE_ISO_LJ(X, POTLIST, N_ATOM_PARTNERS, POTSCALE, PARAMS, ENERGY, GRAD, GTEST, STEST) 
267:  
268:     ! Because we're passing in the entire coordinates array here, we need to know the total number of atoms 
269:     USE COMMONS, ONLY: NATOMS 
270: !    USE MODHESS 
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)  :: POTSCALE             ! The energy unit for this potential, which multiplies E, G and Hess. 
278:     DOUBLE PRECISION, INTENT(IN)  :: PARAMS(10)           ! Maximum number of parameters is hardcoded here 
279:     DOUBLE PRECISION, INTENT(OUT) :: ENERGY               ! The energy of the configuration 
280:     DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOMS)       ! The energy gradient 
281:     LOGICAL, INTENT(IN)           :: GTEST, STEST         ! Flags to specify whether the gradient should be calculated (GTEST) and 
282:                                                           ! whether the Hessian should be calculated (STEST) 
283:  
284:     ! Various powers of the distance between the atoms, and the atom radius 
285:     DOUBLE PRECISION :: DIST, R2(NATOMS,NATOMS), R6, R8(NATOMS,NATOMS), R14(NATOMS,NATOMS), SIG, SIG6, SIG12 
286:     DOUBLE PRECISION :: G(NATOMS,NATOMS), F(NATOMS,NATOMS)  ! G tensor and F tensor (see ISOTROPIC_GRAD and ISOTROPIC_HESSIAN, below) 
287:     DOUBLE PRECISION :: TMP_ENERGY 
288:     INTEGER :: J1, J2, J3, J4, J5, J6 
289:  
290:     SIG = PARAMS(1) 
291:     SIG6 = SIG**6 
292:     SIG12 = SIG6**2 
293:  
294:     TMP_ENERGY=0.0D0 
295:  
296: !    F(:,:) = 0.0D0 
297:                  
298:     ! 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 
299:     ! of DO loops for each call. 
300:     IF (GTEST) THEN 
301:         IF (STEST) THEN  ! Gradient + Hessian 
302:  
303:             ! Rather than looping over all pairs of atoms, we loop over all interacting pairs as specified by POTLIST. 
304:             ! The first entry in each element of POTLIST contains the index of a principal atom. All other entries contain 
305:             ! the interaction partners for this principal atom. 
306:             DO J6=1,UBOUND(POTLIST,1) 
307:                 J1=POTLIST(J6,1)  ! Get the actual atomic index for this atom 
308:                 J3=3*J1 
309:                 R2(J1,J1)=0.0D0 
310:                 R8(J1,J1)=0.0D0 
311:                 R14(J1,J1)=0.0D0 
312:                 G(J1,J1)=0.0D0 
313:                 F(J1,J1)=0.0D0 
314:  
315:                 ! Loop over the partners for atom J1 
316:                 DO J5=2,N_ATOM_PARTNERS(J6)+1  ! Have to add the +1 because the first element in POTLIST(J1,:) 
317:                                                ! is the atom index J1, instead of a neighbour index. 
318:  
319:                     J2=POTLIST(J6,J5)  ! Get the actual atomic index for this atom 
320:                     ! For historical reasons, the order of the variables J2, J4, J5 is odd. 
321:                     J4=3*J2 
322:                     
323:                     R2(J2,J1)=(X(J3-2)-X(J4-2))**2+(X(J3-1)-X(J4-1))**2+(X(J3)-X(J4))**2 
324:                     R2(J2,J1)=1.0D0/R2(J2,J1) 
325:                     R6=R2(J2,J1)**3 
326:  
327:                     TMP_ENERGY=TMP_ENERGY+SIG6*R6*(SIG6*R6-1.0D0) 
328:  
329: !                    write(*,*) "E+G+H. Pair, Dist, Energy:", J1, J2, 1.0D0/R2(J2,J1), SIG6*R6*(SIG6*R6-1.0D0) + 0.25D0 
330:                     ! Set up storage arrays to use in the gradient- and hessian-calculating routines 
331:                     R8(J2,J1)=R2(J2,J1)**4 
332:                     R14(J2,J1)=R8(J2,J1)*R8(J2,J1)/R2(J2,J1) 
333:                     R2(J1,J2)=R2(J2,J1) 
334:                     G(J2,J1)=-24.0D0*(2.0D0*SIG6*R6-1.0D0)*R2(J1,J2)*SIG6*R6 
335:                     G(J1,J2)=G(J2,J1) 
336:                     F(J2,J1)=672.0D0*R14(J2,J1)*SIG12-192.0D0*R8(J2,J1)*SIG6 
337:                     F(J1,J2)=F(J2,J1) 
338:                 ENDDO 
339:             ENDDO 
340:         ELSE             ! Energy + Gradient only 
341:             DO J6=1,UBOUND(POTLIST,1) 
342:                 J1=POTLIST(J6,1) 
343:                 J3=3*J1 
344:                 G(J1,J1)=0.0D0 
345:  
346:                 DO J5=2,N_ATOM_PARTNERS(J6)+1 
347:                     
348:                     J2=POTLIST(J6,J5)  ! For historical reasons, the order of the variables J2, J4, J5 is odd. 
349:                     J4=3*J2 
350:  
351:                     DIST=(X(J3-2)-X(J4-2))**2+(X(J3-1)-X(J4-1))**2+(X(J3)-X(J4))**2 
352:                     DIST=1.0D0/DIST 
353:                     R6=DIST**3 
354:  
355:                     TMP_ENERGY=TMP_ENERGY+R6*SIG6*(R6*SIG6-1.0D0) 
356: !                    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 
357:                     G(J2,J1)=-24.0D0*(2.0D0*R6*SIG6-1.0D0)*DIST*R6*SIG6 
358:                     G(J1,J2)=G(J2,J1) 
359:                 ENDDO 
360:             ENDDO 
361:         ENDIF 
362:     ELSE                ! Energy only 
363:         DO J6=1,UBOUND(POTLIST,1) 
364:             J1=POTLIST(J6,1) 
365:             J3=3*J1 
366: !            G(J1,J1)=0.0D0 
367:  
368:             DO J5=2,N_ATOM_PARTNERS(J6)+1 
369:  
370:                 J2=POTLIST(J6,J5)  ! For historical reasons, the order of the variables J2, J4, J5 is odd. 
371:                 J4=3*J2 
372:  
373:                 DIST=(X(J3-2)-X(J4-2))**2+(X(J3-1)-X(J4-1))**2+(X(J3)-X(J4))**2 
374:                 DIST=1.0D0/DIST 
375:                 R6=DIST**3 
376:  
377:                 TMP_ENERGY=TMP_ENERGY+R6*SIG6*(R6*SIG6-1.0D0) 
378: !                write(*,*) "E. Pair, Dist, Energy:", J1, J2, R2(J2,J1), SIG6*R6*(SIG6*R6-1.0D0) + 0.25D0 
379:             ENDDO 
380:         ENDDO 
381:  
382:     ENDIF 
383:     TMP_ENERGY=4.0D0*TMP_ENERGY 
384:     ENERGY = ENERGY + POTSCALE*TMP_ENERGY 
385:  
386:     IF (.NOT.GTEST) RETURN 
387:  
388:     ! G should already be set 
389:     CALL EXCLUDE_ISOTROPIC_GRAD(POTLIST, N_ATOM_PARTNERS, X, G, GRAD, POTSCALE) 
390:  
391:     IF (.NOT.STEST) RETURN  ! It is assumed we will never need the Hessian without also needing the gradient. 
392:  
393:     ! This is currently a very inefficient implementation! I haven't implemented the use of POTLISTS yet. 
394:     CALL EXCLUDE_ISOTROPIC_HESSIAN(POTLIST, N_ATOM_PARTNERS, X, G, F, R2, POTSCALE) 
395:  
396:     RETURN 
397:  
398: END SUBROUTINE EXCLUDE_ISO_LJ 
399:  
400: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
401: ! WCA potential with support for excluded interactions. 
402: ! Atom radius sigma is given as a variable (SIG), but will often be set to 1. 
403: SUBROUTINE EXCLUDE_ISO_WCA(X, POTLIST, N_ATOM_PARTNERS, POTSCALE, PARAMS, ENERGY, GRAD, GTEST, STEST) 
404:  
405:     ! Because we're passing in the entire coordinates array here, we need to know the total number of atoms 
406:     USE COMMONS, ONLY: NATOMS 
407: !    USE MODHESS 
408:  
409:     IMPLICIT NONE 
410:  
411:     DOUBLE PRECISION, INTENT(IN)  :: X(3*NATOMS)          ! The full atomistic coordinate array 
412:     INTEGER, INTENT(IN)           :: POTLIST(:,:)         ! An array containing the interaction partners for each atom 
413:     INTEGER, INTENT(IN)           :: N_ATOM_PARTNERS(:)   ! An array containing the number of interaction partners for each atom 
414:     DOUBLE PRECISION, INTENT(IN)  :: POTSCALE             ! The energy unit for this potential, which multiplies E, G and Hess. 
415:     DOUBLE PRECISION, INTENT(IN)  :: PARAMS(10)           ! Maximum number of parameters is hardcoded here 
416:     DOUBLE PRECISION, INTENT(OUT) :: ENERGY               ! The energy of the configuration 
417:     DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOMS)       ! The energy gradient 
418:     LOGICAL, INTENT(IN)           :: GTEST, STEST         ! Flags to specify whether the gradient should be calculated (GTEST) and 
419:                                                           ! whether the Hessian should be calculated (STEST) 
420:  
421:     ! Various powers of the distance between the atoms, and the atom radius 
422:     DOUBLE PRECISION :: DIST, R2(NATOMS,NATOMS), R6, R8(NATOMS,NATOMS), R14(NATOMS,NATOMS), SIG, SIG6, SIG12 
423:     DOUBLE PRECISION :: G(NATOMS,NATOMS), F(NATOMS,NATOMS)  ! G tensor and F tensor (see ISOTROPIC_GRAD and ISOTROPIC_HESSIAN, below) 
424:     DOUBLE PRECISION :: TMP_ENERGY 
425:     INTEGER :: J1, J2, J3, J4, J5, J6 
426:  
427:     SIG = PARAMS(1) 
428:     SIG6 = SIG**6 
429:     SIG12 = SIG6**2 
430:  
431:     TMP_ENERGY=0.0D0 
432: !    G(:,:)=0.0D0 
433: !    F(:,:)=0.0D0 
434:  
435:     ! 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 
436:     ! of DO loops for each call. 
437:     IF (GTEST) THEN 
438:         IF (STEST) THEN  ! Gradient + Hessian 
439:  
440:             ! Rather than looping over all pairs of atoms, we loop over all interacting pairs as specified by POTLIST. 
441:             ! The first entry in each element of POTLIST contains the index of a principal atom. All other entries contain 
442:             ! the interaction partners for this principal atom. 
443:             DO J6=1,UBOUND(POTLIST,1) 
444:                 J1=POTLIST(J6,1)  ! Get the actual atomic index for this atom 
445:  
446:                 J3=3*J1 
447:                 R2(J1,J1)=0.0D0 
448:                 R8(J1,J1)=0.0D0 
449:                 R14(J1,J1)=0.0D0 
450:                 G(J1,J1)=0.0D0 
451:                 F(J1,J1)=0.0D0 
452:  
453:                 ! Loop over the partners for atom J1 
454:                 DO J5=2,N_ATOM_PARTNERS(J6)+1  ! Have to add the +1 because the first element in POTLIST(J1,:) 
455:                                                ! is the atom index J1, instead of a neighbour index. 
456:  
457:                     J2=POTLIST(J6,J5)  ! Get the actual atomic index for this atom 
458:  
459:                     ! For historical reasons, the order of the variables J2, J4, J5 is odd. 
460:                     J4=3*J2 
461:  
462:                     R2(J2,J1)=(X(J3-2)-X(J4-2))**2+(X(J3-1)-X(J4-1))**2+(X(J3)-X(J4))**2 
463:  
464:                     IF(R2(J2,J1).GT.WCA_CUT*SIG*SIG) THEN   ! WCA_CUT is a PARAMETER - see top of file. 
465:                     ! We don't compute the energy for this pair of atoms. Set G and F to 0 so that the gradient and 
466:                     ! Hessian terms will go to 0 also. 
467: !                        write(*,*) "E+G+H. Separation>cutoff:", J1, J2, R2(J2,J1) 
468:                         G(J1,J2) = 0.0D0 
469:                         G(J2,J1) = 0.0D0 
470:                         F(J1,J2) = 0.0D0 
471:                         F(J2,J1) = 0.0D0 
472:                         CYCLE 
473:                     ENDIF 
474:  
475:                     R2(J2,J1)=1.0D0/R2(J2,J1) 
476:                     R6=R2(J2,J1)**3 
477:  
478:                     TMP_ENERGY=TMP_ENERGY+SIG6*R6*(SIG6*R6-1.0D0) + 0.25D0 
479: !                    write(*,*) "E+G+H. Pair, Dist, Energy:", J1, J2, 1.0D0/R2(J2,J1), SIG6*R6*(SIG6*R6-1.0D0) + 0.25D0 
480:                     ! Set up storage arrays to use in the gradient- and hessian-calculating routines 
481:                     R8(J2,J1)=R2(J2,J1)**4 
482:                     R14(J2,J1)=R8(J2,J1)*R8(J2,J1)/R2(J2,J1) 
483:                     R2(J1,J2)=R2(J2,J1) 
484:                     G(J2,J1)=-24.0D0*(2.0D0*SIG6*R6-1.0D0)*R2(J1,J2)*SIG6*R6 
485:                     G(J1,J2)=G(J2,J1) 
486:                     F(J2,J1)=672.0D0*R14(J2,J1)*SIG12-192.0D0*R8(J2,J1)*SIG6 
487:                     F(J1,J2)=F(J2,J1) 
488:                 ENDDO 
489:             ENDDO 
490:         ELSE             ! Energy + Gradient only 
491:             DO J6=1,UBOUND(POTLIST,1) 
492:                 J1=POTLIST(J6,1) 
493:                 J3=3*J1 
494:                 G(J1,J1)=0.0D0 
495:  
496:                 DO J5=2,N_ATOM_PARTNERS(J6)+1 
497:                     
498:                     J2=POTLIST(J6,J5)  ! For historical reasons, the order of the variables J2, J4, J5 is odd. 
499:                     J4=3*J2 
500:  
501:                     DIST=(X(J3-2)-X(J4-2))**2+(X(J3-1)-X(J4-1))**2+(X(J3)-X(J4))**2 
502:  
503:                     IF(DIST.GT.WCA_CUT*SIG*SIG) THEN   ! WCA_CUT is a PARAMETER - see top of file. 
504: !                        write(*,*) "E+G. Separation>cutoff:", J1, J2, DIST 
505:                         G(J1,J2) = 0.0D0 
506:                         G(J2,J1) = 0.0D0 
507:                         CYCLE 
508:                     ENDIF 
509:  
510:                     DIST=1.0D0/DIST 
511:                     R6=DIST**3 
512:  
513:                     TMP_ENERGY=TMP_ENERGY+R6*SIG6*(R6*SIG6-1.0D0) + 0.25D0 
514: !                    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 
515:                     G(J2,J1)=-24.0D0*(2.0D0*R6*SIG6-1.0D0)*DIST*R6*SIG6 
516:                     G(J1,J2)=G(J2,J1) 
517:                 ENDDO 
518:             ENDDO 
519:         ENDIF 
520:     ELSE                ! Energy only 
521:         DO J6=1,UBOUND(POTLIST,1) 
522:             J1=POTLIST(J6,1) 
523:             J3=3*J1 
524:             G(J1,J1)=0.0D0 
525:  
526:             DO J5=2,N_ATOM_PARTNERS(J6)+1 
527:  
528:                 J2=POTLIST(J6,J5)  ! For historical reasons, the order of the variables J2, J4, J5 is odd. 
529:                 J4=3*J2 
530:  
531:                 DIST=(X(J3-2)-X(J4-2))**2+(X(J3-1)-X(J4-1))**2+(X(J3)-X(J4))**2 
532:  
533:                 IF(DIST.GT.WCA_CUT*SIG*SIG) THEN   ! WCA_CUT is a PARAMETER - see top of file. 
534: !                    write(*,*) "E. Separation>cutoff:", J1, J2, R2(J2,J1) 
535:                     CYCLE 
536:                 ENDIF 
537:  
538:                 DIST=1.0D0/DIST 
539:                 R6=DIST**3 
540:  
541:                 TMP_ENERGY=TMP_ENERGY+R6*SIG6*(R6*SIG6-1.0D0) + 0.25D0 
542: !                write(*,*) "E. Pair, Dist, Energy:", J1, J2, R2(J2,J1), SIG6*R6*(SIG6*R6-1.0D0) + 0.25D0 
543:             ENDDO 
544:         ENDDO 
545:  
546:     ENDIF 
547:     TMP_ENERGY=4.0D0*TMP_ENERGY 
548:     ENERGY = ENERGY + POTSCALE*TMP_ENERGY 
549:  
550:     IF (.NOT.GTEST) RETURN 
551:  
552:     ! G should already be set 
553:     CALL EXCLUDE_ISOTROPIC_GRAD(POTLIST, N_ATOM_PARTNERS, X, G, GRAD, POTSCALE) 
554:  
555:     IF (.NOT.STEST) RETURN  ! It is assumed we will never need the Hessian without also needing the gradient. 
556:  
557:     ! This is currently a very inefficient implementation! I haven't implemented the use of POTLISTS yet. 
558:     CALL EXCLUDE_ISOTROPIC_HESSIAN(POTLIST, N_ATOM_PARTNERS, X, G, F, R2, POTSCALE) 
559:  
560:     RETURN 
561:  
562: END SUBROUTINE EXCLUDE_ISO_WCA 
563:  
564: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
565: ! END OF ISOTROPIC POTENTIALS261: ! END OF ISOTROPIC POTENTIALS
566: ! Next, two functions that are general to all isotropic potentials. To be clear, by "isotropic", I mean "depends only262: ! Next, two functions that are general to all isotropic potentials. To be clear, by "isotropic", I mean "depends only
567: ! on the distance between the two particles"263: ! on the distance between the two particles"
568: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!264: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
569: 265: 
570: ! For all isotropic potentials the gradient is calculated in the same way. We simply need to know the positions of the266: ! For all isotropic potentials the gradient is calculated in the same way. We simply need to know the positions of the
571: ! two particles, and the G tensor: G = (1/r)dU/dr for an isotropic potential U(r)267: ! two particles, and the G tensor: G = (1/r)dU/dr for an isotropic potential U(r)
572: SUBROUTINE ISOTROPIC_GRAD(N, X, G, TMP_G)268: SUBROUTINE ISOTROPIC_GRAD(N, X, G, TMP_G)
573:     IMPLICIT NONE269:     IMPLICIT NONE
574:     INTEGER, INTENT(IN)           :: N  ! Number of atoms270:     INTEGER, INTENT(IN)           :: N  ! Number of atoms
579: 275: 
580:     DO J1=1,N  ! The atom for which we are calculating the gradient276:     DO J1=1,N  ! The atom for which we are calculating the gradient
581:         J3=3*J1277:         J3=3*J1
582:         DUMMYX = 0.0D0278:         DUMMYX = 0.0D0
583:         DUMMYY = 0.0D0279:         DUMMYY = 0.0D0
584:         DUMMYZ = 0.0D0280:         DUMMYZ = 0.0D0
585:         DO J4=1,N  ! Consider the interaction with each other atom in turn.281:         DO J4=1,N  ! Consider the interaction with each other atom in turn.
586:             ! This inner loop is the slow bit. Minimise the number of operations required here.282:             ! This inner loop is the slow bit. Minimise the number of operations required here.
587:             J2=3*J4283:             J2=3*J4
588:             XMUL=G(J4,J1)  ! Only do the array-access once to save time284:             XMUL=G(J4,J1)  ! Only do the array-access once to save time
589:  
590:             DUMMYX=DUMMYX+XMUL*(X(J3-2)-X(J2-2))285:             DUMMYX=DUMMYX+XMUL*(X(J3-2)-X(J2-2))
591:             DUMMYY=DUMMYY+XMUL*(X(J3-1)-X(J2-1))286:             DUMMYY=DUMMYY+XMUL*(X(J3-1)-X(J2-1))
592:             DUMMYZ=DUMMYZ+XMUL*(X(J3)-X(J2))287:             DUMMYZ=DUMMYZ+XMUL*(X(J3)-X(J2))
593:         ENDDO288:         ENDDO
594:         TMP_G(J3-2) = DUMMYX289:         TMP_G(J3-2) = DUMMYX
595:         TMP_G(J3-1) = DUMMYY290:         TMP_G(J3-1) = DUMMYY
596:         TMP_G(J3) = DUMMYZ291:         TMP_G(J3) = DUMMYZ
597:     ENDDO292:     ENDDO
598: 293: 
599:     RETURN294:     RETURN
639:         ENDDO334:         ENDDO
640:     ENDDO335:     ENDDO
641: 336: 
642: !  Case III, different atoms, same cartesian coordinate.337: !  Case III, different atoms, same cartesian coordinate.
643: 338: 
644:     DO J1=1,N339:     DO J1=1,N
645:         DO J2=1,3340:         DO J2=1,3
646:             J3=3*(J1-1)+J2341:             J3=3*(J1-1)+J2
647:             DO J4=J1+1,N342:             DO J4=J1+1,N
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)343:                 TMP_HESS(3*(J4-1)+J2,J3)=-F(J4,J1)*R2(J4,J1)*(X(J3)-X(3*(J4-1)+J2))**2-G(J4,J1)
649:             ENDDO344:            ENDDO
650:         ENDDO345:         ENDDO
651:     ENDDO346:     ENDDO
652: 347: 
653: !  Case IV: different atoms and different cartesian coordinates.348: !  Case IV: different atoms and different cartesian coordinates.
654: 349: 
655:     DO J1=1,N350:     DO J1=1,N
656:         DO J2=1,3351:         DO J2=1,3
657:             J3=3*(J1-1)+J2352:             J3=3*(J1-1)+J2
658:             DO J4=J1+1,N353:             DO J4=J1+1,N
659:                 DO J5=1,J2-1354:                 DO J5=1,J2-1
669: 364: 
670:     DO J1=1,3*N365:     DO J1=1,3*N
671:         DO J2=J1+1,3*N366:         DO J2=J1+1,3*N
672:             TMP_HESS(J1,J2)=TMP_HESS(J2,J1)367:             TMP_HESS(J1,J2)=TMP_HESS(J2,J1)
673:         ENDDO368:         ENDDO
674:     ENDDO369:     ENDDO
675: 370: 
676:     RETURN371:     RETURN
677: END SUBROUTINE ISOTROPIC_HESSIAN372: END SUBROUTINE ISOTROPIC_HESSIAN
678: 373: 
679:  
680: ! For all isotropic potentials the gradient is calculated in the same way. We simply need to know the positions of the 
681: ! two particles, and the G tensor: G = (1/r)dU/dr for an isotropic potential U(r) 
682: SUBROUTINE EXCLUDE_ISOTROPIC_GRAD(POTLIST, N_ATOM_PARTNERS, X, G, GRAD, POTSCALE) 
683:  
684:     IMPLICIT NONE 
685:     INTEGER, INTENT(IN)           :: POTLIST(:,:) 
686:     INTEGER, INTENT(IN)           :: N_ATOM_PARTNERS(:)  ! Number of atoms 
687:     DOUBLE PRECISION, INTENT(IN)  :: X(:), G(:,:) 
688:     DOUBLE PRECISION, INTENT(OUT) :: GRAD(:)  ! This is the real full gradient, so we don't want to zero it! 
689:     DOUBLE PRECISION, INTENT(IN)  :: POTSCALE 
690:     INTEGER :: J1, J2, J3, J4, J5, J6 
691:     DOUBLE PRECISION :: DUMMYX, DUMMYY, DUMMYZ, XMUL 
692:  
693:     DO J6=1,UBOUND(POTLIST,1)  ! The atom for which we are calculating the gradient 
694:         J1=POTLIST(J6,1) 
695:         J3=3*J1 
696:  
697:         DUMMYX = 0.0D0 
698:         DUMMYY = 0.0D0 
699:         DUMMYZ = 0.0D0 
700:  
701:         ! This inner loop is the slow bit. Minimise the number of operations required here. 
702:         DO J5=2,N_ATOM_PARTNERS(J6)+1 
703:             J4=POTLIST(J6,J5) 
704:             J2=3*J4 
705:             XMUL=G(J4,J1)  ! Only do the array-access once to save time 
706:  
707:             DUMMYX=XMUL*(X(J3-2)-X(J2-2))*POTSCALE 
708:             DUMMYY=XMUL*(X(J3-1)-X(J2-1))*POTSCALE 
709:             DUMMYZ=XMUL*(X(J3)-X(J2))*POTSCALE 
710:  
711:             GRAD(J3-2) = GRAD(J3-2) + DUMMYX 
712:             GRAD(J3-1) = GRAD(J3-1) + DUMMYY 
713:             GRAD(J3) = GRAD(J3) + DUMMYZ 
714:             ! Symmetrise. The G tensor is symmetric (because U(r) depends only on r and r is symmetric wrt exchange of particles) 
715:             ! So the only thing which changes is X1-X2 becomes X2-X1, i.e. a sign change. 
716:             GRAD(J2-2) = GRAD(J2-2) - DUMMYX 
717:             GRAD(J2-1) = GRAD(J2-1) - DUMMYY 
718:             GRAD(J2) = GRAD(J2) - DUMMYZ 
719:  
720:         ENDDO 
721:     ENDDO 
722:  
723:     RETURN 
724: END SUBROUTINE EXCLUDE_ISOTROPIC_GRAD 
725:  
726: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
727: ! For all isotropic potentials the Hessian is calculated in the same way. We simply need to know the positions of the 
728: ! two particles, the G tensor (see above) and the F tensor: F = r*d[(1/r)dU/dr]/dr 
729: SUBROUTINE EXCLUDE_ISOTROPIC_HESSIAN(POTLIST, N_ATOM_PARTNERS, X, G, F, R2, POTSCALE) 
730: ! NOTE! THIS ROUTINE DOES NOT WORK CURRENTLY! USE NUMERICAL HESSIANS INSTEAD! 
731:  
732:  
733: !    USE COMMONS, ONLY: NATOMS 
734:     USE MODHESS 
735:  
736:     IMPLICIT NONE 
737:  
738:     INTEGER, INTENT(IN)           :: POTLIST(:,:) 
739:     INTEGER, INTENT(IN)           :: N_ATOM_PARTNERS(:)  ! Number of atoms 
740:     DOUBLE PRECISION, INTENT(IN)  :: X(:), G(:,:), F(:,:), R2(:,:), POTSCALE 
741:     INTEGER :: J1, J2, J3, J4, J5, J6, J7, J8 
742:     DOUBLE PRECISION :: DUMMY 
743:  
744:     write(*,*) "Hessian for exclusion potentials not currently working. Use numerical hessians instead." 
745:     STOP 
746:  
747: !       Compute the hessian. First are the entirely diagonal terms (3N of them) 
748: !       These correspond to 2nd derivatives wrt the same coordinate twice 
749:  
750:     DO J6=1,UBOUND(POTLIST,1) 
751:         J1=POTLIST(J6,1) 
752:  
753:         DO J2=1,3 
754:             J3=3*(J1-1)+J2 
755:             DUMMY=0.0D0 
756:  
757:             DO J5=2,N_ATOM_PARTNERS(J6)+1 
758:                 J4=POTLIST(J6,J5) 
759:  
760:                 DUMMY=DUMMY+F(J4,J1)*R2(J4,J1)*(X(J3)-X(3*(J4-1)+J2))**2 + G(J4,J1) 
761:             ENDDO 
762:             HESS(J3,J3)=HESS(J3,J3)+DUMMY*POTSCALE 
763:         ENDDO 
764:     ENDDO 
765:  
766: !  Next are the terms where x_i and x_j are on the same atom 
767: !  but are different, e.g. y and z. 
768:  
769:  
770:     DO J6=1,UBOUND(POTLIST,1) 
771:         J1=POTLIST(J6,1) 
772:  
773:         DO J2=1,3 
774:             J3=3*(J1-1)+J2 
775:             DO J4=J2+1,3 
776:                 DUMMY=0.0D0 
777:  
778:                 DO J7=2,N_ATOM_PARTNERS(J6)+1 
779:                     J5=POTLIST(J6,J7) 
780:  
781:                     DUMMY=DUMMY + F(J5,J1)*R2(J5,J1)*(X(J3)-X(3*(J5-1)+J2))*(X(3*(J1-1)+J4)-X(3*(J5-1)+J4)) 
782:                 ENDDO 
783:                 HESS(3*(J1-1)+J4,J3)=HESS(3*(J1-1)+J4,J3)+DUMMY*POTSCALE 
784:                 HESS(J3,3*(J1-1)+J4)=HESS(J3,3*(J1-1)+J4)+DUMMY*POTSCALE ! Symmetrise 
785:             ENDDO 
786:         ENDDO 
787:     ENDDO 
788:  
789: !  Case III, different atoms, same cartesian coordinate. 
790:  
791:  
792:     DO J6=1,UBOUND(POTLIST,1) 
793:         J1=POTLIST(J6,1) 
794:  
795:         DO J2=1,3 
796:             J3=3*(J1-1)+J2 
797:             DUMMY=0.0D0 
798:  
799:             DO J5=2,N_ATOM_PARTNERS(J6)+1 
800:                 J4=POTLIST(J6,J5) 
801:  
802:                 DUMMY=-F(J4,J1)*R2(J4,J1)*(X(J3)-X(3*(J4-1)+J2))**2-G(J4,J1) 
803:                 HESS(3*(J4-1)+J2,J3)=HESS(3*(J4-1)+J2,J3)+DUMMY*POTSCALE 
804:                 HESS(J3,3*(J4-1)+J2)=HESS(J3,3*(J4-1)+J2)+DUMMY*POTSCALE  ! Symmetrise 
805:             ENDDO            
806:         ENDDO 
807:     ENDDO 
808:  
809: !  Case IV: different atoms and different cartesian coordinates. 
810:  
811:     DO J7=1,UBOUND(POTLIST,1) 
812:         J1=POTLIST(J7,1) 
813:  
814:         DO J2=1,3 
815:             J3=3*(J1-1)+J2 
816:  
817:             DO J8=2,N_ATOM_PARTNERS(J7)+1 
818:                 J4=POTLIST(J7,J8) 
819:  
820:                 DO J5=1,J2-1 
821:                     J6=3*(J4-1)+J5 
822:                     DUMMY = -F(J4,J1)*R2(J4,J1)*(X(J3)-X(3*(J4-1)+J2))*(X(3*(J1-1)+J5)-X(J6)) 
823:                     HESS(J6,J3)=HESS(J6,J3)+DUMMY*POTSCALE 
824:                     HESS(J3,J6)=HESS(J3,J6)+DUMMY*POTSCALE ! Symmetrise 
825:                     HESS(3*(J4-1)+J2,3*(J1-1)+J5)=HESS(3*(J4-1)+J2,3*(J1-1)+J5)+DUMMY*POTSCALE 
826:                     HESS(3*(J1-1)+J5,3*(J4-1)+J2)=HESS(3*(J1-1)+J5,3*(J4-1)+J2)+DUMMY*POTSCALE ! Symmetrise 
827:                 ENDDO 
828:             ENDDO 
829:         ENDDO 
830:     ENDDO 
831:  
832:     RETURN 
833: END SUBROUTINE EXCLUDE_ISOTROPIC_HESSIAN 
834:  
835:  
836: END MODULE ISOTROPIC_POTENTIALS374: END MODULE ISOTROPIC_POTENTIALS


r30623/keywords.f 2016-07-06 15:37:52.709078910 +0100 r30622/keywords.f 2016-07-06 15:37:54.933108971 +0100
4614:          MSTRANST=.TRUE.4614:          MSTRANST=.TRUE.
4615: 4615: 
4616:       ELSE IF (WORD.EQ.'MULLERBROWN') THEN4616:       ELSE IF (WORD.EQ.'MULLERBROWN') THEN
4617:          MULLERBROWNT=.TRUE.4617:          MULLERBROWNT=.TRUE.
4618: 4618: 
4619:       ELSE IF (WORD.EQ.'MULTIPLICITY') THEN4619:       ELSE IF (WORD.EQ.'MULTIPLICITY') THEN
4620:          CALL READI(XMUL)4620:          CALL READI(XMUL)
4621: 4621: 
4622:       ELSE IF (WORD .EQ. 'MULTIPOT') THEN4622:       ELSE IF (WORD .EQ. 'MULTIPOT') THEN
4623:          MULTIPOTT = .TRUE.4623:          MULTIPOTT = .TRUE.
 4624:          CALL MULTIPOT_INITIALISE        
4624:          4625:          
4625: !4626: !
4626: ! NOT DOCUMENTED4627: ! NOT DOCUMENTED
4627: !4628: !
4628:       ELSE IF (WORD.EQ.'MYSD') THEN4629:       ELSE IF (WORD.EQ.'MYSD') THEN
4629:          MYSDT=.TRUE.4630:          MYSDT=.TRUE.
4630:          CALL READF(SDTOL)4631:          CALL READF(SDTOL)
4631: 4632: 
4632:       ELSE IF (WORD.EQ.'NATB') THEN4633:       ELSE IF (WORD.EQ.'NATB') THEN
4633:          NATBT=.TRUE.4634:          NATBT=.TRUE.


r30623/main.F 2016-07-06 15:37:53.053083561 +0100 r30622/main.F 2016-07-06 15:37:55.269113510 +0100
 25:       USE PERMU 25:       USE PERMU
 26:       USE F1COM 26:       USE F1COM
 27:       USE MODAMBER 27:       USE MODAMBER
 28:       USE MODAMBER9, only : AMBFINALIO_NODE,MDCRD_UNIT,MDINFO_UNIT,AMBPDB_UNIT, ATMASS1 28:       USE MODAMBER9, only : AMBFINALIO_NODE,MDCRD_UNIT,MDINFO_UNIT,AMBPDB_UNIT, ATMASS1
 29:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_MASSES 29:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_MASSES
 30:       USE MODCHARMM 30:       USE MODCHARMM
 31:       USE PORFUNCS 31:       USE PORFUNCS
 32:       USE TWIST_MOD 32:       USE TWIST_MOD
 33:       USE HOMOREFMOD 33:       USE HOMOREFMOD
 34:       USE GENRIGID, only : RIGIDINIT, GENRIGID_READ_FROM_FILE, DEGFREEDOMS 34:       USE GENRIGID, only : RIGIDINIT, GENRIGID_READ_FROM_FILE, DEGFREEDOMS
 35:       USE MULTIPOT, only: MULTIPOT_INITIALISE 
 36:       USE GAUSS_MOD, ONLY: KEGEN 35:       USE GAUSS_MOD, ONLY: KEGEN
 37:       IMPLICIT NONE 36:       IMPLICIT NONE
 38:       !EXTERNAL READ_CMD_ARGS 37:       !EXTERNAL READ_CMD_ARGS
 39: #ifdef MPI 38: #ifdef MPI
 40:       INCLUDE 'mpif.h' 39:       INCLUDE 'mpif.h'
 41: #endif 40: #endif
 42:       INTEGER J1,J2, JP, MPIERR, NDUMMY3,NPTOTAL,VERSIONTEMP,GETUNIT,LUNIT 41:       INTEGER J1,J2, JP, MPIERR, NDUMMY3,NPTOTAL,VERSIONTEMP,GETUNIT,LUNIT
 43:       DOUBLE PRECISION, ALLOCATABLE :: SCREENC(:) 42:       DOUBLE PRECISION, ALLOCATABLE :: SCREENC(:)
 44:       DOUBLE PRECISION POTEL, DUMMY, INTFREEZETOLSAVE, RMAT(3,3), DUMMY2, DIST2, LINTCONSTRAINTTOL 43:       DOUBLE PRECISION POTEL, DUMMY, INTFREEZETOLSAVE, RMAT(3,3), DUMMY2, DIST2, LINTCONSTRAINTTOL
 45:       DOUBLE PRECISION, ALLOCATABLE :: TMPCOORDS(:), ENDCOORDS(:,:) 44:       DOUBLE PRECISION, ALLOCATABLE :: TMPCOORDS(:), ENDCOORDS(:,:)
188:           CALL GENRIGID_READ_FROM_FILE ()187:           CALL GENRIGID_READ_FROM_FILE ()
189: ! csw34> Tell the user how many degrees of freedom there are in the system188: ! csw34> Tell the user how many degrees of freedom there are in the system
190:           WRITE(MYUNIT, '(A,I15)') " genrigid> rbodyconfig used to specifiy rigid bodies, degrees of freedom now ", DEGFREEDOMS189:           WRITE(MYUNIT, '(A,I15)') " genrigid> rbodyconfig used to specifiy rigid bodies, degrees of freedom now ", DEGFREEDOMS
191:           IF (GCBHT) THEN190:           IF (GCBHT) THEN
192: ! csw34> Make sure we aren't running GCBH with rigid bodies - very dangerous! Could be done but only with great care...191: ! csw34> Make sure we aren't running GCBH with rigid bodies - very dangerous! Could be done but only with great care...
193:               WRITE(MYUNIT, '(A)') " genrigid> ERROR: cannot use rigid bodies with GCBH! Stopping."192:               WRITE(MYUNIT, '(A)') " genrigid> ERROR: cannot use rigid bodies with GCBH! Stopping."
194:               STOP                193:               STOP                
195:           END IF194:           END IF
196:       END IF195:       END IF
197: 196: 
198:       IF (MULTIPOTT) THEN 
199:           CALL MULTIPOT_INITIALISE() 
200:       ENDIF 
201:  
202:       IF (CUDAT) THEN197:       IF (CUDAT) THEN
203:          IF ((CUDAPOT .EQ. 'A') .AND. (.NOT. AMBER12T)) THEN198:          IF ((CUDAPOT .EQ. 'A') .AND. (.NOT. AMBER12T)) THEN
204:             WRITE(MYUNIT,'(A)') " main> The AMBER12 keyword must be used with 'CUDA A'. "199:             WRITE(MYUNIT,'(A)') " main> The AMBER12 keyword must be used with 'CUDA A'. "
205:             STOP200:             STOP
206:          END IF201:          END IF
207:          IF (DEBUG) THEN202:          IF (DEBUG) THEN
208:             CUDAFILENAME1 ="GPU_debug_out"203:             CUDAFILENAME1 ="GPU_debug_out"
209:             OPEN(UNIT=325, FILE=CUDAFILENAME1, STATUS="REPLACE")204:             OPEN(UNIT=325, FILE=CUDAFILENAME1, STATUS="REPLACE")
210:             CLOSE(325)205:             CLOSE(325)
211:          END IF206:          END IF


r30623/multipot.f90 2016-07-06 15:37:54.241099539 +0100 r30622/multipot.f90 2016-07-06 15:37:56.273126940 +0100
  2:   2: 
  3: ! Checklist for adding a new potential:  3: ! Checklist for adding a new potential:
  4: !   a) Write a subroutine to calculate the energy (preferably saving it in isotropic_potentials.f90 or pairwise_potentials.f90)  4: !   a) Write a subroutine to calculate the energy (preferably saving it in isotropic_potentials.f90 or pairwise_potentials.f90)
  5: !   b) Make sure the signature of your subroutine matches the interface in COMPUTE_PAIRWISE_POTENTIAL or COMPUTE_ISOTROPIC_POTENTIAL  5: !   b) Make sure the signature of your subroutine matches the interface in COMPUTE_PAIRWISE_POTENTIAL or COMPUTE_ISOTROPIC_POTENTIAL
  6: !      (see below). Is neither of these is appropriate for your potential, you will need to write a new COMPUTE subroutine.  6: !      (see below). Is neither of these is appropriate for your potential, you will need to write a new COMPUTE subroutine.
  7: !   c) Add your potential to the SELECT CASE statement in MULTIPOT_CALL  7: !   c) Add your potential to the SELECT CASE statement in MULTIPOT_CALL
  8: !   d) If your potential requires atom indices to be input in any format other than a simple list, add your requirements to the  8: !   d) If your potential requires atom indices to be input in any format other than a simple list, add your requirements to the
  9: !      SELECT CASE statement in MULTIPOT_INITIALISE  9: !      SELECT CASE statement in MULTIPOT_INITIALISE
 10: !   e) Check that you've made equivalent changes in GMIN so that the two systems stay in sync. 10: !   e) Check that you've made equivalent changes in GMIN so that the two systems stay in sync.
 11: !   f) Enjoy! 11: !   f) Enjoy!
  12: 
 12: MODULE MULTIPOT 13: MODULE MULTIPOT
 13:  14: 
 14:     USE PAIRWISE_POTENTIALS 15:     USE PAIRWISE_POTENTIALS
 15:     USE ISOTROPIC_POTENTIALS 16:     USE ISOTROPIC_POTENTIALS
 16:     USE COMMONS, ONLY: DEBUG 17:     USE COMMONS, ONLY: DEBUG
 17:  18: 
 18:     IMPLICIT NONE 19:     IMPLICIT NONE
 19:  20: 
 20:     ! We pre-calculate as much as possible, hence the need for lots of different storage arrays, documented here. 21:     ! We pre-calculate as much as possible, hence the need for lots of different storage arrays, documented here.
 21:     INTEGER :: NPOTTYPES = 0                        ! The number of different potential types being used 22:     INTEGER :: NPOTTYPES = 0                        ! The number of different potential types being used
 33:     INTEGER, ALLOCATABLE :: POTLISTS(:,:,:)         ! Holds the indices of all atoms belonging to each potential type, and 34:     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 35:                                                     ! 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 36:                                                     ! 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 37:                                                     ! type POTTYPES(m), and POTLISTS(m,n,2:N_ATOM_PARTNERS(m,n)) are the indices of
 37:                                                     ! its partners. 38:                                                     ! its partners.
 38:  39: 
 39:     INTEGER, ALLOCATABLE :: DEGS_BY_POT(:,:)        ! For use by isotropic potentials, where all atoms interact with all other atoms. 40:     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 41:                                                     ! Contains a list of the degrees of freedom belonging to all the atoms using this
 41:                                                     ! potential. 42:                                                     ! potential.
 42:  43: 
 43:  
 44: CONTAINS 44: CONTAINS
 45: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 45: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 46:     ! 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
 47: SUBROUTINE MULTIPOT_INITIALISE 47: SUBROUTINE MULTIPOT_INITIALISE
 48:  48: 
 49:     USE COMMONS, ONLY: NATOMS 49:     USE COMMONS, ONLY: NATOMS
 50:     USE GENRIGID, ONLY: RIGIDINIT, RB_BY_ATOM 50:     USE GENRIGID, ONLY: RIGIDINIT, RB_BY_ATOM
 51:  51: 
 52:     IMPLICIT NONE 52:     IMPLICIT NONE
 53:  53: 
 54:     INTEGER NPARAMS, ATOM1, ATOM2 54:     INTEGER NPARAMS, ATOM1, ATOM2
 55:     INTEGER :: ATOMLIST(NATOMS) 55:     INTEGER :: ATOMLIST(NATOMS)
 56:     CHARACTER(LEN=1000) :: DUMMYCHAR 56:     CHARACTER(LEN=1000) :: DUMMYCHAR
 57:     LOGICAL :: END 57:     LOGICAL :: END
 58:     INTEGER :: J1, J2, J3, J4, iostatus, COUNTER 58:     INTEGER :: J1, J2, J3, iostatus, COUNTER
 59:  
 60:     ! Variables needed to read in exclusion lists 
 61:     INTEGER :: N_EXCLUDE_LINES, MAX_LINE_LENGTH 
 62:     INTEGER, ALLOCATABLE :: LINE_LEN(:) 
 63:     INTEGER, ALLOCATABLE :: EXCLUSIONS(:,:)  ! Holds the indices of excluded interactions. EXCLUSIONS(l,m) contains 
 64:                                              ! the m'th atom in the l'th set of excluded interactions for the current potential. 
 65:  59: 
 66:     ! Input file: multipotconfig 60:     ! Input file: multipotconfig
 67:     ! Format as follows. 61:     ! Format as follows.
 68:     ! 62:     !
 69:     ! Introduce each type of interaction with a line: 'POT' 63:     ! Introduce each type of interaction with a line: 'POT'
 70:     ! To comment out a potential type, simply insert a '#' at the start of the 'POT' line. 64:     ! To comment out a potential type, simply insert a '#' at the start of the 'POT' line.
 71:     ! The next line has the form: POTTYPE n scale nparams 65:     ! The next line has the form: POTTYPE n scale nparams
 72:     !   POTTYPE is a string specifying the type of potential 66:     !   POTTYPE is a string specifying the type of potential
 73:     !   n is the number of atoms which will use this type of potential 67:     !   n is the number of atoms which will use this type of potential
 74:     !   scale is the factor by which the energy of this interaction must be scaled to match the others 68:     !   scale is the factor by which the energy of this interaction must be scaled to match the others
 75:     !   nparams is the number of arguments to the pairwise potential, which will be potential specific. nparams<=10, currently. 69:     !   nparams is the number of arguments to the pairwise potential, which will be potential specific. nparams<=10, currently.
 76:     ! The next line is the list of input parameters (up to 10 of them), separated by spaces.  70:     ! The next line is the list of input parameters (up to 10 of them), separated by spaces. 
 77:     !   Note that nparams needs to be at least 1. If your potential has no extra parameters, set nparams to 1 71:     !   Note that nparams needs to be at least 1. If your potential has no extra parameters, set nparams to 1
 78:     !   and have this line contain a single "0", which you then don't need to use. 72:     !   and have this line contain a single "0", which you then don't need to use.
 79:     ! The next line(s) contain lists of the atoms which will interact according to the current potential. 73:     ! The next line(s) contain lists of the atoms which will interact according to the current potential.
 80:     !   See below for example formats which may be used. 74:     !   See below for example formats which may be used.
 81:     ! Some formats (notably EWCA, ELJ) may require additional input lines in multipotconfig. 
 82:  75: 
 83: !   Read multipotconfig once to count the number of types of potential 76: !   Read multipotconfig once to count the number of types of potential
 84:     NPOTTYPES=0 77:     NPOTTYPES=0
 85:     OPEN(UNIT=22,FILE='multipotconfig',status='old') 78:     OPEN(UNIT=22,FILE='multipotconfig',status='old')
 86:     DO 79:     DO
 87:        READ(22,*,IOSTAT=iostatus) DUMMYCHAR 80:        READ(22,*,IOSTAT=iostatus) DUMMYCHAR
 88:        IF (iostatus<0) THEN 81:        IF (iostatus<0) THEN
 89:           CLOSE(22) 82:           CLOSE(22)
 90:           EXIT 83:           EXIT
 91:        ELSE IF (DUMMYCHAR(1:3).EQ.'POT') THEN 84:        ELSE IF (DUMMYCHAR(1:3).EQ.'POT') THEN
132:         IF (DUMMYCHAR(1:1).EQ.'#') THEN  ! This potential is commented out: skip it.125:         IF (DUMMYCHAR(1:1).EQ.'#') THEN  ! This potential is commented out: skip it.
133:             IF(DEBUG) WRITE(*,*) "Skipping commented block in multipotconfig"126:             IF(DEBUG) WRITE(*,*) "Skipping commented block in multipotconfig"
134:             DO127:             DO
135:                 READ(22,*) DUMMYCHAR128:                 READ(22,*) DUMMYCHAR
136:                 IF (DUMMYCHAR(1:3).EQ.'POT') EXIT  ! We have found the start of the next uncommented header129:                 IF (DUMMYCHAR(1:3).EQ.'POT') EXIT  ! We have found the start of the next uncommented header
137:             ENDDO130:             ENDDO
138:         ENDIF131:         ENDIF
139: 132: 
140:         READ(22,*) DUMMYCHAR133:         READ(22,*) DUMMYCHAR
141:         READ(22,*) DUMMYCHAR134:         READ(22,*) DUMMYCHAR
142:  
143:  
144:         ! Now read in the atom numbers from the following line(s).135:         ! Now read in the atom numbers from the following line(s).
145:         ! 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: 
146:         SELECT CASE(POTTYPES(J1))138:         SELECT CASE(POTTYPES(J1))
147: 139: 
148:         ! 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.
149:         ! This method of reading input should extend simply to 3- and 4-body potentials (e.g. dihedral terms).141:         CASE('HSPR')
150:         CASE('HSPR', 'PWCA', 'PLJ')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)
151:             DO J2 = 1, NATOM_BY_POT(J1)/2  ! J2 ranges over the number of pairs (Note, NATOM_BY_POT contains the actual no of atoms)143:                 ! We only want to compute the potential once for each spring.
152:                 ! We only want to compute the potential once for each pair. 
153: 144: 
154:                 READ(22,*) ATOM1, ATOM2145:                 READ(22,*) ATOM1, ATOM2
155: 146: 
156:                 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
157:                 POTLISTS(J1,2*J2,1) = ATOM2148:                 POTLISTS(J1,2*J2,1) = ATOM2
158: 149: 
159:                 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...
160:                 POTLISTS(J1,2*J2-1,2) = ATOM2151:                 POTLISTS(J1,2*J2-1,2) = ATOM2
161:                 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)
162:             ENDDO153:             ENDDO
163: 154: 
164:             IF(DEBUG) THEN155:             IF(DEBUG) THEN
165:                 WRITE(*,*) "Potential:", POTTYPES(J1)156:                 WRITE(*,*) "Potential:", POTTYPES(J1)
166:                 WRITE(*,*) "N_ATOM_PARTNERS:"157:                 WRITE(*,*) "Bonds:"
167:                 WRITE(*,*) N_ATOM_PARTNERS(J1,:NATOM_BY_POT(J1)) 
168:                 WRITE(*,*) "POTLISTS:" 
169:                 DO J2 = 1,NATOM_BY_POT(J1)158:                 DO J2 = 1,NATOM_BY_POT(J1)
170:                     WRITE(*,*) "Atom number", J2, "in this potential"159:                     IF (N_ATOM_PARTNERS(J1,J2).EQ.1) WRITE(*,*) POTLISTS(J1,J2,1:2)
171:                     WRITE(*,*) POTLISTS(J1,J2,:N_ATOM_PARTNERS(J1,J2)) 
172:                 ENDDO160:                 ENDDO
173:             ENDIF161:             ENDIF
174: 162: 
175:         ! For "iso_" potentials. Every specified atom with this potential type interacts with all the others.163:         ! For "iso_" potentials. Every specified atom interacts with all the others, but this is done by means of a single
176:         ! The input format is simple: list all the atoms using this potential on a single line, separated by spaces.164:         ! call to a hard-coded version of the potential (e.g. ISO_LJ) rather than a series of calls to PAIRWISE_LJ for every pair of
 165:         ! interacting atoms. This is less flexible, but significantly more efficient for large systems.
 166:         ! The input format is straightforward: simply list all the atoms using this potential on a single line, separated by spaces.
177:         ! They don't have to be in index order, but everything will make more sense if they are!167:         ! They don't have to be in index order, but everything will make more sense if they are!
178:         ! There is currently no facility to exclude the interactions within a rigid body. Use EWCA/ELJ if you need that. 
179:         CASE('ILJ','IWCA')168:         CASE('ILJ','IWCA')
180:             READ(22,*) ATOMLIST(:NATOM_BY_POT(J1))169:             READ(22,*) ATOMLIST(:NATOM_BY_POT(J1))
181: 170: 
182:             ! All we need to save for this type of potential is the list of whole-system degrees of freedom on which171:             ! All we need to save for this type of potential is the list of whole-system degrees of freedom on which
183:             ! the potential will depend.172:             ! the potential will depend.
184:             DO J2=1,NATOM_BY_POT(J1)173:             DO J2=1,NATOM_BY_POT(J1)
185:                 J3 = ATOMLIST(J2)174:                 J3 = ATOMLIST(J2)
186:                 DEGS_BY_POT(J1,3*J2-2) = 3*J3-2175:                 DEGS_BY_POT(J1,3*J2-2) = 3*J3-2
187:                 DEGS_BY_POT(J1,3*J2-1) = 3*J3-1176:                 DEGS_BY_POT(J1,3*J2-1) = 3*J3-1
188:                 DEGS_BY_POT(J1,3*J2)   = 3*J3177:                 DEGS_BY_POT(J1,3*J2)   = 3*J3
189:             ENDDO178:             ENDDO
190: 179: 
191:             IF(DEBUG) THEN180:             IF(DEBUG) THEN
192:                 write(*,*) "Potential:", POTTYPES(J1)181:                 write(*,*) "Potential:", POTTYPES(J1)
193:                 write(*,*) "DEGS_BY_POT:"182:                 write(*,*) "DEGS_BY_POT:"
194:                 write(*,*) DEGS_BY_POT(J1,:NATOM_BY_POT(J1)*3)183:                 write(*,*) DEGS_BY_POT(J1,:NATOM_BY_POT(J1)*3)
195:             ENDIF184:             ENDIF
196: 185: 
197:         CASE('EWCA', 'ELJ')186:         CASE DEFAULT
198: 187:         ! Default: an isotropic potential handled pair-by-pair. NOTE: likely to be very inefficient for large systems. Use the "iso_"
199:         ! An isotropic potential with the option to switch off particular interactions.188:         ! version of the relevant potential instead (which has the same input format).
200:         !189:         ! All atoms with this potential type interact with all others - unless they are in the same rigid body.
201:         ! Specify the allowed interactions as normal for isotropic potentials: with a list of atom indices on a single line.190:         ! All atoms interacting according to that potential should be specified on a single line, separated by spaces
202:         ! If RIGIDINIT is set, then interactions within a rigid body are not calculated.191:         ! They don't have to be in index order, but everything will make more sense if they are!
203:         ! Additional excluded interactions are specified by appending the following to this multipotconfig entry: 
204:         ! 
205:         ! EXCLUDE nlines max_length 
206:         ! line_len_1 
207:         ! index1 index2 ... 
208:         ! line_len_2 
209:         ! index1 index2 ... 
210:         ! ... 
211:         ! 
212:         ! nlines is the number of sets of interactions you wish to exclude 
213:         ! max_length is the largest number of atoms involved in an excluded set (interactions will not be calculated for any 
214:         ! pair of atoms within the set) 
215:         ! The next nlines pairs of lines specify each excluded set. line_len_1 specifies the number of atoms in excluded set 1, 
216:         ! and the following line contains the indices of these atoms. 
217:         ! 
218:         ! If this potential type is the last one in your MULTIPOTCONFIG file, then you must include a final line  
219:  
220:             ! Read the list of interacting atoms as normal 
221:             READ(22,*) ATOMLIST(:NATOM_BY_POT(J1))192:             READ(22,*) ATOMLIST(:NATOM_BY_POT(J1))
222: 193: 
223:             ! Parse the exclusions lists. Read the header line once to check that it is an EXCLUDE header line as expected.194:             IF(RIGIDINIT .AND. .FALSE.) THEN   ! We can avoid calculate interactions within a rigid body by only including atom partner pairs
224:             READ(22,*) DUMMYCHAR195:                                                ! that involve atoms from two different bodies
225:             IF (DUMMYCHAR(1:7) .EQ. 'EXCLUDE') THEN 
226:                 ! Re-read this line to extract the parameters that we want 
227:                 BACKSPACE(22) 
228:                 READ(22,*) DUMMYCHAR, N_EXCLUDE_LINES, MAX_LINE_LENGTH 
229:  
230:                 ALLOCATE(LINE_LEN(N_EXCLUDE_LINES), EXCLUSIONS(N_EXCLUDE_LINES,MAX_LINE_LENGTH)) 
231:                             
232:                 ! Read in the lists of atoms for which interactions will be excluded 
233:                     DO J2 = 1, N_EXCLUDE_LINES 
234:                     READ(22,*) LINE_LEN(J2)  ! This line tells us how many indices to read from the next line 
235:                     READ(22,*) EXCLUSIONS(J2,:LINE_LEN(J2)) 
236:                 ENDDO 
237:  
238:             ELSE 
239:                 WRITE(*,*) "multipot> Warning. An exclusion potential is being used with no EXCLUDE lists specified." 
240:                 WRITE(*,*) "multipot> It's probably more efficient to use the ISO- version of this potential instead." 
241:                 BACKSPACE(22)  ! Step back one record so that we can move on to the next potential 
242:             ENDIF 
243:  
244:             ! We do not calculate interactions which have been explicitly excluded, or interactions within rigid bodies. 
245:             ! If there are no rigid bodies and no exclusions have been specified, we can use a simpler method to construct 
246:             ! POTLISTS (see the ELSE block) 
247:             IF(RIGIDINIT .OR. ALLOCATED(EXCLUSIONS)) THEN   
248:  
249:                 ! We need the list of which rigid body each atom belongs to. GENRIGID_READ_FROM_FILE must already have been 
250:                 ! called. For GMIN, this IF block should never be satisfied but for OPTIM the RIGIDINIT keyword needs to be 
251:                 ! placed higher in the odata file than MULTIPOT. 
252:                 IF (.NOT. ALLOCATED(RB_BY_ATOM)) THEN 
253:                     WRITE(*,*) "multipot> Error: RB_BY_ATOM is not allocated by the time MULTIPOT_INITIALISE is called." 
254:                     WRITE(*,*) "multipot> The RIGIDINIT keyword must precede the MULTIPOT keyword in your odata file." 
255:                     STOP 98 
256:                 ENDIF 
257:                                  
258:                 ! Work out which interactions are still allowed after accounting for exclusions, and put them into POTLISTS. 
259:                 DO J2=1,NATOM_BY_POT(J1)196:                 DO J2=1,NATOM_BY_POT(J1)
260:                     POTLISTS(J1,J2,1) = ATOMLIST(J2)   ! The first element in POTLISTS(x,y,:) is the index of atom y in the list197:                     POTLISTS(J1,J2,1) = ATOMLIST(J2)   ! Remember, the first element in POTLISTS(x,y,:) is the index of atom y in the list
261:                                                        ! for potential type x. The remaining elements are the partners of atom y.198:                                                        ! for potential type x. The remaining elements are the partners of atom y.
262:                     N_ATOM_PARTNERS(J1,J2) = 0199:                     N_ATOM_PARTNERS(J1,J2) = 0
263: 200:                     DO J3=J2,NATOM_BY_POT(J1)
264:                     ! Work out which other atoms in this potential will interact with POTLISTS(J1,J2,1).201:                     ! So that each pair only gets calculated once, each partner pair is only included once. So the 2nd atom
265:                     ! Only include each pair once to avoid double counting, so the following loop runs over indices higher than J2.202:                     ! in each pair to appear in the list gets listed as a partner to the 1st atom, but not the other way round.
266:                     inner: DO J3=J2+1,NATOM_BY_POT(J1)203:                         IF (RB_BY_ATOM(ATOMLIST(J3)) .NE. RB_BY_ATOM(ATOMLIST(J2))) THEN
267:                     204:                             ! Only add partners that belong to different rigid bodies
268:                         IF ((.NOT. RIGIDINIT) .OR. (RB_BY_ATOM(ATOMLIST(J3)) .NE. RB_BY_ATOM(ATOMLIST(J2)))) THEN 
269:                             ! Only add partners that belong to different rigid bodies  
270:                             ! (always fulfilled if there are no rigid bodies!) 
271:  
272:                             IF(ALLOCATED(EXCLUSIONS)) THEN 
273:                                 ! Check to see if this interaction was excluded. 
274:                                 DO J4=1,N_EXCLUDE_LINES 
275:                                     IF ((ANY(EXCLUSIONS(J4,:)==ATOMLIST(J2))) .AND. (ANY(EXCLUSIONS(J4,:)==ATOMLIST(J3)))) THEN 
276:                                         CYCLE inner  ! Don't incude this interaction. 
277:                                     ENDIF 
278:                                 ENDDO 
279:                             ENDIF 
280:  
281:                             ! If we have got this far, this pair must be a genuine interacting pair. 
282:                             ! Avoid double counting: each interacting pair is only included once.  
283:                             ! (The fact that J3 starts from J2+1 ensures this) 
284:                             N_ATOM_PARTNERS(J1,J2) = N_ATOM_PARTNERS(J1,J2) + 1205:                             N_ATOM_PARTNERS(J1,J2) = N_ATOM_PARTNERS(J1,J2) + 1
285:                             POTLISTS(J1,J2,1+N_ATOM_PARTNERS(J1,J2)) = ATOMLIST(J3)206:                             POTLISTS(J1,J2,1+N_ATOM_PARTNERS(J1,J2)) = ATOMLIST(J3)
286:                         ENDIF207:                         ENDIF
287:                     ENDDO inner208:                     ENDDO
288:                 ENDDO209:                 ENDDO
289:             ELSE   210:             ELSE   ! No rigid bodies so all atoms in ATOMLIST should partner with all others.
290:                ! No rigid bodies or exclusions so all atoms in ATOMLIST should partner with all others. 
291:                 DO J2 = 1,NATOM_BY_POT(J1)  ! We still only add each pair once.211:                 DO J2 = 1,NATOM_BY_POT(J1)  ! We still only add each pair once.
292:                     N_ATOM_PARTNERS(J1,J2) = NATOM_BY_POT(J1)-J2  ! So each atom partners with every atom that comes later in the list212:                     N_ATOM_PARTNERS(J1,J2) = NATOM_BY_POT(J1)-J2  ! So each atom partners with every atom that comes later in the list
293:                     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))
294:                 ENDDO214:                 ENDDO
295:             ENDIF        215:             ENDIF
296: 216: 
297:             IF(DEBUG) THEN217:             IF(DEBUG) THEN
298:                 WRITE(*,*) "Potential:", POTTYPES(J1)218:                 WRITE(*,*) "Potential:", POTTYPES(J1)
299:                 WRITE(*,*) "Number of atoms:", NATOM_BY_POT(J1) 
300:                 WRITE(*,*) "N_ATOM_PARTNERS:"219:                 WRITE(*,*) "N_ATOM_PARTNERS:"
301:                 WRITE(*,*) N_ATOM_PARTNERS(J1,:NATOM_BY_POT(J1))220:                 WRITE(*,*) N_ATOM_PARTNERS(J1,:NATOM_BY_POT(J1))
302:                 WRITE(*,*) "POTLISTS:"221:                 WRITE(*,*) "POTLISTS:"
303:                 DO J2 = 1,NATOM_BY_POT(J1)222:                 DO J2 = 1,NATOM_BY_POT(J1)
304:                     WRITE(*,*) "Atom number", J2, "in this potential"223:                     IF (N_ATOM_PARTNERS(J1,J2).GT.0) THEN
305:                     WRITE(*,*) POTLISTS(J1,J2,:1+N_ATOM_PARTNERS(J1,J2))224:                        WRITE(*,*) "Atom number", J2, "in this potential"
 225:                        WRITE(*,*) POTLISTS(J1,J2,N_ATOM_PARTNERS(J1,J2))
 226:                     ENDIF
306:                 ENDDO227:                 ENDDO
307:             ENDIF228:             ENDIF
308: 229: 
309:         CASE DEFAULT 
310:            WRITE(*,*) "multipot> Error: no rule to read input for this type of potential." 
311:            STOP 
312:         END SELECT230:         END SELECT
313:     ENDDO231:     ENDDO
314:     CLOSE(22)232:     CLOSE(22)
315: 233: 
316: END SUBROUTINE MULTIPOT_INITIALISE234: END SUBROUTINE MULTIPOT_INITIALISE
317:  
318: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!235: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
319: ! Calculate the potential energy and the gradient (if GTEST is set) and hessian (if STEST)236: ! Calculate the potential energy and the gradient (if GTEST is set) and hessian (if STEST)
320: SUBROUTINE MULTIPOT_CALL (X, G, ENERGY, GTEST, STEST)237: SUBROUTINE MULTIPOT_CALL (X, G, ENERGY, GTEST, STEST)
321:     USE COMMONS, ONLY: NATOMS238:     USE COMMONS, ONLY: NATOMS
322:     USE MODHESS239:     USE MODHESS
323: 240: 
324:     IMPLICIT NONE241:     IMPLICIT NONE
325: 242: 
326:     DOUBLE PRECISION, INTENT(IN)  :: X(3*NATOMS)243:     DOUBLE PRECISION, INTENT(IN)  :: X(3*NATOMS)
327:     DOUBLE PRECISION, INTENT(OUT) :: ENERGY, G(3*NATOMS)244:     DOUBLE PRECISION, INTENT(OUT) :: ENERGY, G(3*NATOMS)
328:     LOGICAL, INTENT(IN)           :: GTEST, STEST245:     LOGICAL, INTENT(IN)           :: GTEST, STEST
329:     INTEGER                       :: J1246:     INTEGER                       :: J1
330: 247: 
331:     INTEGER                       :: THISN 
332:     DOUBLE PRECISION              :: THESEPARAMS(10) 
333:  
334:     ENERGY  = 0.D0248:     ENERGY  = 0.D0
335:     IF (GTEST) G(:) = 0.D0249:     IF (GTEST) G(:) = 0.D0
336:     IF (STEST) THEN250:     IF (STEST) THEN
337:         HESS(:,:) = 0.D0251:         HESS(:,:) = 0.D0
338:     ENDIF252:     ENDIF
339: 253: 
340:     ! Cycle through the different types of interaction and perform the potential call for all the corresponding atoms254:     ! Cycle through the different types of interaction and perform the potential call for all the corresponding atoms
341:     DO J1 = 1, NPOTTYPES255:     DO J1 = 1, NPOTTYPES
342:  
343:         ! Extract a few sub-arrays manually to avoid ifort's irritating "array temporary" warning message. 
344:         THISN = NATOM_BY_POT(J1) 
345:         THESEPARAMS = POTPARAMS(J1,:) 
346:  
347:         SELECT CASE(POTTYPES(J1))256:         SELECT CASE(POTTYPES(J1))
348:             CASE('ILJ')257:             CASE('ILJ')
349:                 ! Only one parameter for this potential: the particle radius, sigma258:                 ! Only one parameter for this potential: the particle radius, sigma
350:                 CALL COMPUTE_ISOTROPIC_POTENTIAL(X, G, ENERGY, GTEST, STEST, ISO_LJ, J1, THISN)259:                 CALL COMPUTE_ISOTROPIC_POTENTIAL(X, G, ENERGY, GTEST, STEST, ISO_LJ, J1, NATOM_BY_POT(J1))
351:             CASE('IWCA')260:             CASE('IWCA')
352:                 ! Only one parameter for this potential: the particle radius, sigma261:                 ! Only one parameter for this potential: the particle radius, sigma
353:                 CALL COMPUTE_ISOTROPIC_POTENTIAL(X, G, ENERGY, GTEST, STEST, ISO_WCA, J1, THISN)262:                 CALL COMPUTE_ISOTROPIC_POTENTIAL(X, G, ENERGY, GTEST, STEST, ISO_WCA, J1, NATOM_BY_POT(J1))
354:             CASE('PLJ')263:             CASE('PLJ')
355:                 ! Only one parameter for this potential: the particle radius, sigma264:                 ! Only one parameter for this potential: the particle radius, sigma
356:                 CALL COMPUTE_PAIRWISE_POTENTIAL(X, G, ENERGY, GTEST, STEST, PAIRWISE_LJ, J1)265:                 CALL COMPUTE_PAIRWISE_POTENTIAL(X, G, ENERGY, GTEST, STEST, PAIRWISE_LJ, J1)
357:             CASE('PWCA')266:             CASE('PWCA')
358:                 ! Only one parameter for this potential: the particle radius, sigma267:                 ! Only one parameter for this potential: the particle radius, sigma
359:                 CALL COMPUTE_PAIRWISE_POTENTIAL(X, G, ENERGY, GTEST, STEST, PAIRWISE_WCA, J1)268:                 CALL COMPUTE_PAIRWISE_POTENTIAL(X, G, ENERGY, GTEST, STEST, PAIRWISE_WCA, J1)
360:             CASE('HSPR')269:             CASE('HSPR')
361:                 ! Parameter is the equilibrium bond length, R0. Energy is returned in units of k_spr.270:                 ! Parameter is the equilibrium bond length, R0. Energy is returned in units of k_spr.
362:                 CALL COMPUTE_PAIRWISE_POTENTIAL(X, G, ENERGY, GTEST, STEST, HARMONIC_SPRINGS, J1)271:                 CALL COMPUTE_PAIRWISE_POTENTIAL(X, G, ENERGY, GTEST, STEST, HARMONIC_SPRINGS, J1)
363:             ! For exclusion potentials, we must also pass in a list specifying the number of interacting partners possessed 
364:             ! by each atom, and a list specifying which atoms make up these partners (from POTLISTS) 
365:             CASE('ELJ') 
366:                 ! Only one parameter for this potential: the particle radius, sigma 
367:                 CALL EXCLUDE_ISO_LJ(X, POTLISTS(J1,:THISN,:), N_ATOM_PARTNERS(J1,:THISN), POTSCALES(J1), THESEPARAMS, ENERGY, G, GTEST, STEST) 
368:             CASE('EWCA') 
369:                  ! Only one parameter for this potential: the particle radius, sigma 
370:                 CALL EXCLUDE_ISO_WCA(X, POTLISTS(J1,:THISN,:), N_ATOM_PARTNERS(J1,:THISN), POTSCALES(J1), THESEPARAMS, ENERGY, G, GTEST, STEST) 
371:             CASE DEFAULT272:             CASE DEFAULT
372:                 ! We shouldn't every get here, unless you have added a new type of potential to MULTIPOT_INITIALISE and forgotten 
373:                 ! to add it to this SELECT CASE. 
374:                 WRITE(*,*) "multipot> Error: unspecified potential"273:                 WRITE(*,*) "multipot> Error: unspecified potential"
375:                 STOP274:                 STOP
376:         END SELECT275:         END SELECT
377:     ENDDO276:     ENDDO
378: 277: 
 278: 
379: END SUBROUTINE MULTIPOT_CALL279: END SUBROUTINE MULTIPOT_CALL
380: 280: 
381: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!281: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
382: !   Evaluate the energy(+gradient(+hessian)) for a set of atoms interacting according to this particular potential282: !   Evaluate the energy(+gradient(+hessian)) for a set of atoms interacting according to this particular potential
383: SUBROUTINE COMPUTE_ISOTROPIC_POTENTIAL(X, G, ENERGY, GTEST, STEST, POT, POTID, TMP_NATOMS)283: SUBROUTINE COMPUTE_ISOTROPIC_POTENTIAL(X, G, ENERGY, GTEST, STEST, POT, POTID, TMP_NATOMS)
384:     USE COMMONS, ONLY: NATOMS284:     USE COMMONS, ONLY: NATOMS
385:     USE MODHESS285:     USE MODHESS
386:     IMPLICIT NONE286:     IMPLICIT NONE
387:     DOUBLE PRECISION, INTENT(IN)    :: X(3*NATOMS)287:     DOUBLE PRECISION, INTENT(IN)    :: X(3*NATOMS)
388:     DOUBLE PRECISION, INTENT(INOUT) :: G(3*NATOMS)288:     DOUBLE PRECISION, INTENT(INOUT) :: G(3*NATOMS)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0