hdiff output

r32795/commons.f90 2017-06-14 20:30:10.843914868 +0100 r32794/commons.f90 2017-06-14 20:30:13.891954803 +0100
 82:      &                 KINT, INTFREEZETOL, IMSEPMIN, IMSEPMAX, CONCUTABS, CONCUTFRAC, & 82:      &                 KINT, INTFREEZETOL, IMSEPMIN, IMSEPMAX, CONCUTABS, CONCUTFRAC, &
 83:      &                 LPDGEOMDIFFTOL, INTCONFRAC, MAXCONE, INTRMSTOL, BFGSTSTOL, ORBITTOL, & 83:      &                 LPDGEOMDIFFTOL, INTCONFRAC, MAXCONE, INTRMSTOL, BFGSTSTOL, ORBITTOL, &
 84:      &                 INTCONSTRAINTTOL, INTCONSTRAINTDEL, RBCUTOFF, INTCONSTRAINTREP, INTCONSTRAINREPCUT, & 84:      &                 INTCONSTRAINTTOL, INTCONSTRAINTDEL, RBCUTOFF, INTCONSTRAINTREP, INTCONSTRAINREPCUT, &
 85:      &                 INTLJTOL, INTLJDEL, INTLJEPS, REPCON, INTDGUESS, CHECKREPCUTOFF, INTMINFAC, FREEZETOL, & 85:      &                 INTLJTOL, INTLJDEL, INTLJEPS, REPCON, INTDGUESS, CHECKREPCUTOFF, INTMINFAC, FREEZETOL, &
 86:      &                 LOCALPERMCUT, LOCALPERMCUT2, INTCONCUT, QCIRADSHIFT, MLPLAMBDA, & 86:      &                 LOCALPERMCUT, LOCALPERMCUT2, INTCONCUT, QCIRADSHIFT, MLPLAMBDA, &
 87:      &                 CAPSIDRHO,CAPSIDEPS,SIGMAPENT,RADPENT,SIGMAHEX,RADHEX,SIGMAPH, KLIM, SCA, & 87:      &                 CAPSIDRHO,CAPSIDEPS,SIGMAPENT,RADPENT,SIGMAHEX,RADHEX,SIGMAPH, KLIM, SCA, &
 88:      &                 QCIADDREPCUT, QCIADDREPEPS, MLQLAMBDA, TANHFAC, LJADDCUTOFF,LJADDREFNORM 88:      &                 QCIADDREPCUT, QCIADDREPEPS, MLQLAMBDA, TANHFAC, LJADDCUTOFF,LJADDREFNORM
 89:  89: 
 90:       LOGICAL DEBUG, TARGET, MORSET, CUTT, SEEDT, CENT, TSALLIST, FREEZECORE, NEWJUMP, RENORM, CAPSID, FREEZE, & 90:       LOGICAL DEBUG, TARGET, MORSET, CUTT, SEEDT, CENT, TSALLIST, FREEZECORE, NEWJUMP, RENORM, CAPSID, FREEZE, &
 91:      &        OTPT, LJMFT, STRANDT, PAHT, SWT, MSTRANST, STOCKT, STICKYT, BLNT, MYSDT, FREEZERES, CENTXY, & 91:      &        OTPT, LJMFT, STRANDT, PAHT, SWT, MSTRANST, STOCKT, STICKYT, BLNT, MYSDT, FREEZERES, CENTXY, &
 92:      &        MSORIGT, SQUEEZET, PERIODIC, SCT, MSCT, MGUPTAT, RESIZET, TIP, RIGID, CALCQT, MPIT, GBHT, JMT, LJCOULT, LJ_GAUSST, OPPT, SETCENT, & 92:      &        MSORIGT, SQUEEZET, PERIODIC, SCT, MSCT, MGUPTAT, RESIZET, TIP, RIGID, CALCQT, MPIT, GBHT, JMT, LJCOULT, LJ_GAUSST, SETCENT, &
 93:      &        SORTT, HIT, SAVEQ, PARALLELT, FIXD, RKMIN, BSMIN, PERMDIST, PERMOPT, BSWL, BSPT, BSPTRESTART, & 93:      &        SORTT, HIT, SAVEQ, PARALLELT, FIXD, RKMIN, BSMIN, PERMDIST, PERMOPT, BSWL, BSPT, BSPTRESTART, &
 94:      &        SYMMETRIZE, SYMMETRIZECSM, PRINT_PTGRP, PRINT_MINDATA, DUMPT, NEON, ARGON, P46, NORESET, TABOOT, EVSTEPT, PACHECO, DL_POLY, QUCENTRE, & 94:      &        SYMMETRIZE, SYMMETRIZECSM, PRINT_PTGRP, PRINT_MINDATA, DUMPT, NEON, ARGON, P46, NORESET, TABOOT, EVSTEPT, PACHECO, DL_POLY, QUCENTRE, &
 95:      &        STAR, PLUS, TWOPLUS, GROUND, DIPOLE, DFTBT, DFTBCT, SW, SUPERSTEP, EAMLJT, PBGLUET, TRACKDATAT, & 95:      &        STAR, PLUS, TWOPLUS, GROUND, DIPOLE, DFTBT, DFTBCT, SW, SUPERSTEP, EAMLJT, PBGLUET, TRACKDATAT, &
 96:      &        EAMALT, ALGLUET, MGGLUET, GUPTAT, LJATT, FST, DECAY, COOP, FIXBIN, GAUSST, QUENCHDOS, FIXDIHEFLAG, & 96:      &        EAMALT, ALGLUET, MGGLUET, GUPTAT, LJATT, FST, DECAY, COOP, FIXBIN, GAUSST, QUENCHDOS, FIXDIHEFLAG, &
 97:      &        FRAUSIT, ANGST, SELFT, STEPOUT, WENZEL, THRESHOLDT, THOMSONT, MULLERBROWNT, CHARMMENERGIES, & 97:      &        FRAUSIT, ANGST, SELFT, STEPOUT, WENZEL, THRESHOLDT, THOMSONT, MULLERBROWNT, CHARMMENERGIES, &
 98:      &        PROJ, RGCL2, TOSI, WELCH, AXTELL, AMBER, FIXIMAGE, BINARY, SHIFTCUT, ARNO, TUNNELT, TWOD, & 98:      &        PROJ, RGCL2, TOSI, WELCH, AXTELL, AMBER, FIXIMAGE, BINARY, SHIFTCUT, ARNO, TUNNELT, TWOD, &
 99:      &        BLJCLUSTER, BLJCLUSTER_NOCUT, COMPRESST, FIX, FIXT, BFGS, LBFGST, DBRENTT, DZTEST, FNI, FAL, CPMD, TNT, ZETT1, & 99:      &        BLJCLUSTER, BLJCLUSTER_NOCUT, COMPRESST, FIX, FIXT, BFGS, LBFGST, DBRENTT, DZTEST, FNI, FAL, CPMD, TNT, ZETT1, &
100:      &        ZETT2, GBH_RESTART, RESTART, CONJG, NEWRESTART, AVOID, NATBT, DIFFRACTT, CHRMMT, INTMINT, LB2T, &100:      &        ZETT2, GBH_RESTART, RESTART, CONJG, NEWRESTART, AVOID, NATBT, DIFFRACTT, CHRMMT, INTMINT, LB2T, &
101:      &        PTMC, BINSTRUCTURES, PROGRESS, MODEL1T, NEWRESTART_MD, CHANGE_TEMP, NOCISTRANS, CHECKCHIRALITY, &101:      &        PTMC, BINSTRUCTURES, PROGRESS, MODEL1T, NEWRESTART_MD, CHANGE_TEMP, NOCISTRANS, CHECKCHIRALITY, &
102:      &        GBT, GBDT, GBDPT, GEMT, LINRODT, RADIFT, CAPBINT, DBPT, DBPTDT, DMBLMT, DMBLPYT, EFIELDT, PAHAT, STOCKAAT, MORSEDPT, &102:      &        GBT, GBDT, GBDPT, GEMT, LINRODT, RADIFT, CAPBINT, DBPT, DBPTDT, DMBLMT, DMBLPYT, EFIELDT, PAHAT, STOCKAAT, MORSEDPT, &
107:      &        PROJIT, PROJIHT, LEAPDIHE, DUMPQUT, DUMPBESTT, LJSITE, BLJSITE, LJSITEATTR, DUMPSTEPST, PAHAGENRIGIDT, &107:      &        PROJIT, PROJIHT, LEAPDIHE, DUMPQUT, DUMPBESTT, LJSITE, BLJSITE, LJSITEATTR, DUMPSTEPST, PAHAGENRIGIDT, &
108:      &        AMBERT, RANDOMSEEDT, PYGPERIODICT, LJCAPSIDT, PYBINARYT, SWAPMOVEST, MOVABLEATOMST, LIGMOVET,DUMPSTRUCTURES, &108:      &        AMBERT, RANDOMSEEDT, PYGPERIODICT, LJCAPSIDT, PYBINARYT, SWAPMOVEST, MOVABLEATOMST, LIGMOVET,DUMPSTRUCTURES, &
109:      &        LJSITECOORDST, VGW, ACKLANDT, G46, DF1T, PULLT, LOCALSAMPLET, CSMT, A9INTET, INTERESTORE, COLDFUSION, &109:      &        LJSITECOORDST, VGW, ACKLANDT, G46, DF1T, PULLT, LOCALSAMPLET, CSMT, A9INTET, INTERESTORE, COLDFUSION, &
110:      &        CSMGUIDET, MULTISITEPYT, CHAPERONINT, AVOIDRESEEDT, OHCELLT, UNFREEZEFINALQ, PERCOLATET, PERCT, PERCACCEPTED, PERCCOMPMARKOV, PERCGROUPT, &110:      &        CSMGUIDET, MULTISITEPYT, CHAPERONINT, AVOIDRESEEDT, OHCELLT, UNFREEZEFINALQ, PERCOLATET, PERCT, PERCACCEPTED, PERCCOMPMARKOV, PERCGROUPT, &
111:      &        GENALT, MINDENSITYT, RESTRICTREGION, RESTRICTREGIONTEST, RESTRICTCYL, ACK1, ACK2, HARMONICF, PERCGROUPRESEEDT, &111:      &        GENALT, MINDENSITYT, RESTRICTREGION, RESTRICTREGIONTEST, RESTRICTCYL, ACK1, ACK2, HARMONICF, PERCGROUPRESEEDT, &
112:      &        HARMONICDONTMOVE, DUMPUNIQUE, FREEZESAVE, TBP, RBSYMT, PTMCDUMPSTRUCT, PTMCDUMPENERT, PYCOLDFUSION, MONITORT,&112:      &        HARMONICDONTMOVE, DUMPUNIQUE, FREEZESAVE, TBP, RBSYMT, PTMCDUMPSTRUCT, PTMCDUMPENERT, PYCOLDFUSION, MONITORT,&
113:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, POLYT, SANDBOXT, &113:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, POLYT, SANDBOXT, &
114:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &114:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &
115:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &115:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &
116:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &116:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &
117:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, LJ_VAR_RADT, READMASST, SPECMASST, NEWTSALLIST, &117:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &
118:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, &118:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, &
119:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &119:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &
120:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &120:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &
121:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &121:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &
122:      &        MLP3T, MKTRAPT, MLPB3T, MLPB3NEWT, MULTIPOTT, QCIAMBERT, MLPNEWREG, DJWRBT, STEALTHYT, LJADDT, QCINOREPINT, RIGIDMDT, &122:      &        MLP3T, MKTRAPT, MLPB3T, MLPB3NEWT, MULTIPOTT, QCIAMBERT, MLPNEWREG, DJWRBT, STEALTHYT, LJADDT, QCINOREPINT, RIGIDMDT, &
123:      &        DUMPMQT, MLQT, MLQPROB, LJADD2T, MLPVB3T, NOREGBIAS, PYADDT, PYADD2T, LJADD3T, REORDERADDT,  LJADD4T, &123:      &        DUMPMQT, MLQT, MLQPROB, LJADD2T, MLPVB3T, NOREGBIAS, PYADDT, PYADD2T, LJADD3T, REORDERADDT,  LJADD4T, &
124:      &        SQNMT, SQNM_DEBUGT, SQNM_BIOT, BENZRIGIDEWALDT, ORTHO, EWALDT124:      &        SQNMT, SQNM_DEBUGT, SQNM_BIOT, BENZRIGIDEWALDT, ORTHO, EWALDT
125: !125: !
126:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:)126:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:)
127:       DOUBLE PRECISION, ALLOCATABLE :: ATMASS(:)127:       DOUBLE PRECISION, ALLOCATABLE :: ATMASS(:)


r32795/finalio.f90 2017-06-14 20:30:11.063917750 +0100 r32794/finalio.f90 2017-06-14 20:30:14.139958052 +0100
 25:        AMBER12_WRITE_XYZ 25:        AMBER12_WRITE_XYZ
 26:   USE OPEP_INTERFACE_MOD, ONLY: OPEP_FINISH, OPEP_WRITE_PDB 26:   USE OPEP_INTERFACE_MOD, ONLY: OPEP_FINISH, OPEP_WRITE_PDB
 27:   USE QMODULE 27:   USE QMODULE
 28:   USE MODCHARMM 28:   USE MODCHARMM
 29:   USE AMHGLOBALS, ONLY:NMRES,IRES 29:   USE AMHGLOBALS, ONLY:NMRES,IRES
 30:   USE BGUPMOD 30:   USE BGUPMOD
 31:   USE PERMU 31:   USE PERMU
 32:   USE MODHESS, ONLY : HESS, MASSWT 32:   USE MODHESS, ONLY : HESS, MASSWT
 33:   USE CONVEX_POLYHEDRA_MODULE, ONLY: VIEW_POLYHEDRA 33:   USE CONVEX_POLYHEDRA_MODULE, ONLY: VIEW_POLYHEDRA
 34:   USE LJ_GAUSS_MOD, ONLY: VIEW_LJ_GAUSS 34:   USE LJ_GAUSS_MOD, ONLY: VIEW_LJ_GAUSS
 35:   USE OPP_MOD, ONLY: VIEW_OPP 
 36:   USE ORBITALS_MOD, ONLY: ORBITALS_FINISH 35:   USE ORBITALS_MOD, ONLY: ORBITALS_FINISH
 37:  36: 
 38:   IMPLICIT NONE 37:   IMPLICIT NONE
 39:  38: 
 40:   !   MCP 39:   !   MCP
 41:   INTEGER III, I3,  GLY_COUNT, ID, NUMCRD, NUMPRO, NCPHST, GETUNIT, AMHUNIT1, LUNIT, NSYMOPS, NCLOSE 40:   INTEGER III, I3,  GLY_COUNT, ID, NUMCRD, NUMPRO, NCPHST, GETUNIT, AMHUNIT1, LUNIT, NSYMOPS, NCLOSE
 42:   INTEGER J1, J2, J3, J4, J5, MYUNIT2, I1, NDUMMY, MYUNIT3, NRBS1, NRBS2, LJGOUNIT, MJ2, REORDERUNIT 41:   INTEGER J1, J2, J3, J4, J5, MYUNIT2, I1, NDUMMY, MYUNIT3, NRBS1, NRBS2, LJGOUNIT, MJ2, REORDERUNIT
 43:   DOUBLE PRECISION RBCOORDS(NRBSITES*3), DCOORDS(3*NATOMS), EDUMMY, ITDET, DIST 42:   DOUBLE PRECISION RBCOORDS(NRBSITES*3), DCOORDS(3*NATOMS), EDUMMY, ITDET, DIST
 44:   ! 43:   !
 45:   !ds656> Symmetry detection for clusters on a substrate... 44:   !ds656> Symmetry detection for clusters on a substrate...
1612:   ELSE IF (POLYT) THEN1611:   ELSE IF (POLYT) THEN
1613: 1612: 
1614:      CALL VIEW_POLYHEDRA()1613:      CALL VIEW_POLYHEDRA()
1615:      RETURN1614:      RETURN
1616: 1615: 
1617:   ELSE IF (LJ_GAUSST) THEN1616:   ELSE IF (LJ_GAUSST) THEN
1618: 1617: 
1619:      CALL VIEW_LJ_GAUSS()1618:      CALL VIEW_LJ_GAUSS()
1620:      RETURN1619:      RETURN
1621: 1620: 
1622:   ELSE IF (OPPT) THEN 
1623:  
1624:      CALL VIEW_OPP() 
1625:      RETURN 
1626:  
1627:   ELSE IF (PTSTSTT) THEN1621:   ELSE IF (PTSTSTT) THEN
1628: 1622: 
1629:      CALL VIEWPTSTST()1623:      CALL VIEWPTSTST()
1630:      RETURN1624:      RETURN
1631: 1625: 
1632:      !|gd351>1626:      !|gd351>
1633: 1627: 
1634:   ELSE IF (PATCHY) THEN1628:   ELSE IF (PATCHY) THEN
1635: 1629: 
1636:      CALL VIEWPATCHY()1630:      CALL VIEWPATCHY()


r32795/GMIN.tex 2017-06-14 20:30:10.611911826 +0100 r32794/GMIN.tex 2017-06-14 20:30:13.671951920 +0100
575: centre-of-mass distance for each atom.575: centre-of-mass distance for each atom.
576: The compression term can be turned off once a given RMS tolerance is achieved576: The compression term can be turned off once a given RMS tolerance is achieved
577: using the {\it GUIDE} keyword.577: using the {\it GUIDE} keyword.
578: This combination can be used with the generalised rigid body scheme introduced via the578: This combination can be used with the generalised rigid body scheme introduced via the
579: {\it RIGIDINIT} keyword.579: {\it RIGIDINIT} keyword.
580: 580: 
581: \item {\it COMPRESSRIGID kcomp\/ dist}: if at least one rigid body's centre of mass is further than {\it dist} from its nearest neighbour,  581: \item {\it COMPRESSRIGID kcomp\/ dist}: if at least one rigid body's centre of mass is further than {\it dist} from its nearest neighbour,  
582: add a harmonic compression potential with force constant {\it kcomp\/} using the centre-of-mass distance for each rigid body. The compression is 582: add a harmonic compression potential with force constant {\it kcomp\/} using the centre-of-mass distance for each rigid body. The compression is 
583: disabled once the rigid bodies are all within {\it dist} of another. 583: disabled once the rigid bodies are all within {\it dist} of another. 
584: This keyword is only implemented for a few specific potentials. 584: This keyword is only implemented for a few specific potentials. 
585: However, the {\it COMPRESS} and {\it GUIDE} keywords can be used with the general rigid body scheme585: However, the {\it COMPRESS} and {\it GUIDE} keywords can be used with the general ridig body scheme
586: introcduced via {\it GENRIGID}.586: introcduced via {\it GENRIGID}.
587: 587: 
588: \item {\it COMMENT \/}: the rest of the line is ignored.588: \item {\it COMMENT \/}: the rest of the line is ignored.
589: 589: 
590: \item {\it CONCUTABS x \/}: used in the QCI procedure. This is the590: \item {\it CONCUTABS x \/}: used in the QCI procedure. This is the
591: absolute distance constraints are allowed to deviate before the constraint591: absolute distance constraints are allowed to deviate before the constraint
592: spring potential is applied.592: spring potential is applied.
593: Defailt value is 0.15.593: Defailt value is 0.15.
594: 594: 
595: \item {\it COOPMOVE n cut\/}: specifies cooperative moves in the step-taking routine. An atom is595: \item {\it COOPMOVE n cut\/}: specifies cooperative moves in the step-taking routine. An atom is
695: 695: 
696: \item {\it EAMLJ A0 beta Z0\/}: specifies the EAMLJ potential (Baskes, {\it Phys.~Rev.~Lett.\/},696: \item {\it EAMLJ A0 beta Z0\/}: specifies the EAMLJ potential (Baskes, {\it Phys.~Rev.~Lett.\/},
697: {\bf 27}, 2592, 1999) with parameters {\it A0\/}, {\it beta\/} and {\it Z0\/}.697: {\bf 27}, 2592, 1999) with parameters {\it A0\/}, {\it beta\/} and {\it Z0\/}.
698: 698: 
699: \item {\it EDIFF econv\/}: quench minima are only considered to be different if their699: \item {\it EDIFF econv\/}: quench minima are only considered to be different if their
700: energies differ by at least $econv$. This option mainly affects the lowest energy700: energies differ by at least $econv$. This option mainly affects the lowest energy
701: saved geometries. If the current quench energy is within $econv$ of a saved energy, but701: saved geometries. If the current quench energy is within $econv$ of a saved energy, but
702: lies lower, then the saved energy and geometry are replaced.702: lies lower, then the saved energy and geometry are replaced.
703: The default is $0.02$ but different values are appropriate for different potentials.703: The default is $0.02$ but different values are appropriate for different potentials.
704: 704: 
705: \item {\it ENPERMS\/}: Enumerate all distinct permutations of a binary system (0 \textless NTYPEA \textless NATOMS) using Lehmer's algorithm. Instead of a traditional geometry-perturbing step, Lehmer's algorithm will generate a new permutationl isomer by swapping the coordinates of two deterministically-picked unlike atoms. The procedure terminates once all the distinct permutations have been enumerated. 705: \item {\it ENPERMS\/}: Enumerate all distinct permutations of a binary system (0 < NTYPEA < NATOMS) using Lehmer's algorithm. Instead of a traditional geometry-perturbing step, Lehmer's algorithm will generate a new permutationl isomer by swapping the coordinates of two deterministically-picked unlike atoms. The procedure terminates once all the distinct permutations have been enumerated. 
706:  706:  
707: \item {\it EQUILIBRATION equil DumpEveryNthQuench\/}: {\it equil} is the number of 707: \item {\it EQUILIBRATION equil DumpEveryNthQuench\/}: {\it equil} is the number of 
708: MC steps preceding the accumulation of the708: MC steps preceding the accumulation of the
709: density of states histogram in a Wang-Landau709: density of states histogram in a Wang-Landau
710: basin-sampling run. The default is 0. {\it DumpEveryNthQuench} specifies how often the710: basin-sampling run. The default is 0. {\it DumpEveryNthQuench} specifies how often the
711: statistics are recorded into the output files.711: statistics are recorded into the output files.
712: 712: 
713: \item {\it EXEQ n}: for use with the {\it RESERVOIR\/} keyword. After a configuration is713: \item {\it EXEQ n}: for use with the {\it RESERVOIR\/} keyword. After a configuration is
714: taken from the reservoir by the lowest temperature non-reservoir replica we do not714: taken from the reservoir by the lowest temperature non-reservoir replica we do not
715: record visits statistics for {\it n} steps. This allows some local equilibration if715: record visits statistics for {\it n} steps. This allows some local equilibration if
885: {\it GROUP name bondatom1 bondatom2 groupsize rotationscalefactor probselect}885: {\it GROUP name bondatom1 bondatom2 groupsize rotationscalefactor probselect}
886: 886: 
887: {\it groupatom1}887: {\it groupatom1}
888: 888: 
889: {\it groupatom2}889: {\it groupatom2}
890: 890: 
891: {\it groupatom3}891: {\it groupatom3}
892: 892: 
893: \ldots893: \ldots
894: 894: 
895: The group rotation axis is defined by the vector from {\it bondatom1}-\textgreater{\it bondatom2}, and the rotation is scaled by {\it rotationscalefactor}895: The group rotation axis is defined by the vector from {\it bondatom1}->{\it bondatom2}, and the rotation is scaled by {\it rotationscalefactor}
896: .896: .
897: 897: 
898: Here is an example {\textrm atomgroups} file containing two groups:898: Here is an example {\textrm atomgroups} file containing two groups:
899: 899: 
900: {\textrm GROUP OME 6 5 4 1.0 0.8}900: {\textrm GROUP OME 6 5 4 1.0 0.8}
901: 901: 
902: {\textrm 1}902: {\textrm 1}
903: 903: 
904: {\textrm 2}904: {\textrm 2}
905: 905: 
1096: Choosing values according to nearest neighbours in a reference geometry can bias the1096: Choosing values according to nearest neighbours in a reference geometry can bias the
1097: system towards that configuration.1097: system towards that configuration.
1098: 1098: 
1099: \item {\it LJADD2 n\/}: specifies an addressable LJ potential for multiple copies of an $n$-particle1099: \item {\it LJADD2 n\/}: specifies an addressable LJ potential for multiple copies of an $n$-particle
1100: cluster. The data file {\tt epsilon} is required1100: cluster. The data file {\tt epsilon} is required
1101: with $n\times n$ entries, one per line. These well depth parameters are then applied modulo $n$.1101: with $n\times n$ entries, one per line. These well depth parameters are then applied modulo $n$.
1102: For particles with the same address, only the repulsive part of the LJ potential is used.1102: For particles with the same address, only the repulsive part of the LJ potential is used.
1103: 1103: 
1104: 1104: 
1105: \item {\it LJAT $Z^*$ scale\/}: specifies the Lennard-Jones plus Axilrod-Teller 1105: \item {\it LJAT $Z^*$ scale\/}: specifies the Lennard-Jones plus Axilrod-Teller 
1106: potential for a reduced triple-dipole strength of $Z^*$.1106: potential for a reduced triple-dipole struength of $Z^*$.
1107: Provided mostly as a guiding function for {\it DFTBC\/}. 1107: Provided mostly as a guiding function for {\it DFTBC\/}. 
1108: {\it rescale\/} is the coordinate rescaling factor from the LJAT potential1108: {\it rescale\/} is the coordinate rescaling factor from the LJAT potential
1109: to DFTB, default value 2.424.1109: to DFTB, default value 2.424.
1110: 1110: 
1111: \item {\it LJCOUL nc $q'$ f $T_{\rm swap}$} specifies a cluster of Lennard-Jones particles in which the first {\it nc}1111: \item {\it LJCOUL nc $q'$ f $T_{\rm swap}$} specifies a cluster of Lennard-Jones particles in which the first {\it nc}
1112: particles carry identical reduced charges {\it $q'$} in addition to the Lennard-Jones interaction.1112: particles carry identical reduced charges {\it $q'$} in addition to the Lennard-Jones interaction.
1113: The parameter {\it f} specifies what fraction of the Monte Carlo steps should be swaps between the1113: The parameter {\it f} specifies what fraction of the Monte Carlo steps should be swaps between the
1114: positions of a charged and a neutral particle, rather than a conventional step.1114: positions of a charged and a neutral particle, rather than a conventional step.
1115: $T_{\rm swap}$ is the temperature to be used in the acceptance criterion for swap moves, overriding1115: $T_{\rm swap}$ is the temperature to be used in the acceptance criterion for swap moves, overriding
1116: that specified using the {\it TEMPERATURE} keyword.  Generally, a lower temperature is more effective1116: that specified using the {\it TEMPERATURE} keyword.  Generally, a lower temperature is more effective
1117: at finding the lowest-energy permutation of charges.  The default value of $T_{\rm swap}$ is zero.1117: at finding the lowest-energy permutation of charges.  The default value of $T_{\rm swap}$ is zero.
1118: The reduced charge $q'$ is related to the actual charge $q$ by $q'=q/(4\pi\epsilon_0\epsilon\sigma)^{1/2}$,1118: The reduced charge $q'$ is related to the actual charge $q$ by $q'=q/(4\pi\epsilon_0\epsilon\sigma)^{1/2}$,
1119: where $\epsilon$ and $\sigma$ are the Lennard-Jones well depth and length parameter respectively.1119: where $\epsilon$ and $\sigma$ are the Lennard-Jones well depth and length parameter respectively.
1120: This way, the reduced energy of two charges is $E'=q'^2/r'$, where $r'=r/\sigma$ is the reduced distance1120: This way, the reduced energy of two charges is $E'=q'^2/r'$, where $r'=r/\sigma$ is the reduced distance
1121: bdetween the charges.1121: bdetween the charges.
1122: 1122: 
1123: \item {\it LJGAUSS mode rcut $\epsilon$ $r_0$ $\sigma^2$ params} specifies that the Lennard-Jones Gauss double well potential should be used.1123: \item {\it LJGAUSS mode rcut eps $r_0$} sepcifies that the Lennard-Jones Gauss double well potential should be used.
1124: The Lennard-Jones potential has an energy minimum of 1 and a minimum at distance 1.1124: The Lennard-Jones potential has an energy minimum of 1 and a minimum at distance 1.
1125: The relative weight of the Gaussian potential is $\epsilon$ and the Gaussian peak is at distance $r_0$. The width of the Gaussian is1125: The relative weight of the Gaussian potential is {\it eps} and the Gaussian peak is at distance $r_0$. The Stoddard-Ford
1126: controlled by $\sigma^2$. The Stoddard-Ford procedure is used to make the potential and its first derivatives go smoothly to zero at {\it rcut}.1126: procedure is used to make the potential and its first and second derivatives go smoothly to zero at {\it rcut}.
1127: {\it mode} can be either 0, 1, 2 or 3. \textit{mode} 0 specifies the normal potential, and can be either a cluster or with periodic1127: {\it mode} can be either 0, 1 or 2. 0 specifies the normal potential, and can be either a cluster or with periodic
1128: boundary conditions with an orthogonal unit cell using the PERIODIC keyword. \textit{mode} 1 or 2 requires periodic boundary conditions.1128: boundary conditions with an orthogonal unit cell using the PERIODIC keyword. 1 or 2 requires periodic boundary conditions.
1129: The PERIODIC keyword must be used, but the unit cell parameters from it will be ignored. Instead, the last line of the coords file will1129: The PERIODIC keyword must be used, but the unit cell parameters from it will be ignored. Instead, the last line of the coords file will
1130: be interpreted as the unit cell length parameters and will be optimised. Additionally, for {\it mode} 2, the second1130: be interpreted as the unit cell length parameters and will be optimised. Additionally, for {\it mode} 2, the second
1131: last line of the coords file be interpreted as unit cell angle parameters for a triclinic unit cell, and optimised. Since there1131: last line of the coords file be interpreted as unit cell angle parameters for a triclinic unit cell, and optimised. Since there
1132: are combinations of unit cell angles that are physically impossible, steps will be rejected if such combinations are detected.1132: are combinations of unit cell angles that are physically impossible, steps will be rejected if such combinations are detected.
1133: In \textit{mode} 3, as well as the cell being optimised, potential parameters specified by \textit{params} (an integer between 1 and 7) 
1134: will also be optimised. If the 1 bit of \textit{params} is set, $\epsilon$ will be optimised, if the 2 bit is set, $r_0$ will be optimised, 
1135: and if the 4 bit is set, $\sigma^2$ will be optimised. 
1136:  
1137: \item {\it LJVARRAD} specifies that the Lennard-Jones potential be used, where each particle can have a different diameter. The diameter 
1138: of every particle should be specified, one per line, in the file \textit{lj\_var\_rad}, which should be in the working directory. This 
1139: potential is compatible with the generalised rigid-body framework. 
1140: 1133: 
1141: \item {\it LOCALSAMPLE abthresh acthresh\/}: Keyword currently under construction! For three groups of atoms defined in movableatoms1134: \item {\it LOCALSAMPLE abthresh acthresh\/}: Keyword currently under construction! For three groups of atoms defined in movableatoms
1142: (A,B,C), a step is quenched when the centre of coordinates of A-\textgreater B is less than {\it abthresh} AND A-\textgreater C is less than {\it acthresh}. 1135: (A,B,C), a step is quenched when the centre of coordinates of A->B is less than {\it abthresh} AND A->C is less than {\it acthresh}. 
1143: If this condition is broken AFTER the quench, it is automatically rejected.1136: If this condition is broken AFTER the quench, it is automatically rejected.
1144: 1137: 
1145: \item {\it LPERMDIST n thresh cut orbittol\/}: provides the preferred local1138: \item {\it LPERMDIST n thresh cut orbittol\/}: provides the preferred local
1146: permutational alignment procedure.1139: permutational alignment procedure.
1147: {\it n\/} is the maximum number of extra atoms to include in running1140: {\it n\/} is the maximum number of extra atoms to include in running
1148: {\bf minperm\/} for each permutable group (default 4), rather than1141: {\bf minperm\/} for each permutable group (default 4), rather than
1149: treating the whole system at the same time, as for {\it PERMDIST}.1142: treating the whole system at the same time, as for {\it PERMDIST}.
1150: The extra atoms that define the neighbourhood of the permutable group are1143: The extra atoms that define the neighbourhood of the permutable group are
1151: added one at a time in order of their distance from the centre of coordinates of the1144: added one at a time in order of their distance from the centre of coordinates of the
1152: group, averaged over the start and finish geometries.1145: group, averaged over the start and finish geometries.
1403: \item {\it OEINT \/}: interaction energy between 2 peptides will be used as an order parameter---requires 1396: \item {\it OEINT \/}: interaction energy between 2 peptides will be used as an order parameter---requires 
1404: further documentation.1397: further documentation.
1405: 1398: 
1406: \item {\it OHCELL \/}: allow point group operations for a cubic1399: \item {\it OHCELL \/}: allow point group operations for a cubic
1407: supercell in subroutine {\tt minpermdist.f90} permutational distance minimisation.1400: supercell in subroutine {\tt minpermdist.f90} permutational distance minimisation.
1408: 1401: 
1409: \item {\it OPEP option\/}: specifies the use of the OPEP poetential, {\it option\/} is either PROTEIN or RNA. The following additional files are needed:1402: \item {\it OPEP option\/}: specifies the use of the OPEP poetential, {\it option\/} is either PROTEIN or RNA. The following additional files are needed:
1410: {\textrm conf\_initiale.pdb}, which contains the starting configuration and the atom information, {\textrm OPEP\_params}, which can contain additional settings for example for debugging,1403: {\textrm conf\_initiale.pdb}, which contains the starting configuration and the atom information, {\textrm OPEP\_params}, which can contain additional settings for example for debugging,
1411: {\textrm ichain.dat}, {\textrm scale.dat}, {\textrm cutoffs.dat}, {\textrm parametres.top} and {\textrm parametres.list}, which are the generic OPEP files. 1404: {\textrm ichain.dat}, {\textrm scale.dat}, {\textrm cutoffs.dat}, {\textrm parametres.top} and {\textrm parametres.list}, which are the generic OPEP files. 
1412: 1405: 
1413: \item {\it OPP mode rcut k $\phi$ params} specifies that the Oscillating Pair Potential potential should be used. This potential can 
1414: have several wells, depending on the value of \textit{k}, a frequency parameter, and $\phi$, a phase parameter.  
1415: The Stoddard-Ford procedure is used to make the potential and its first derivatives go smoothly to zero at {\it rcut}. 
1416: {\it mode} can be either 0, 1 or 2. \textit{mode} 0 specifies the normal potential, and can be either a cluster or with periodic 
1417: boundary conditions with an orthogonal unit cell using the PERIODIC keyword. \textit{mode} 1 or 2 requires periodic boundary conditions. 
1418: The PERIODIC keyword must be used, but the unit cell parameters from it will be ignored. Instead, the second last line of the coords file 
1419: will be interpreted as unit cell angles and the last line as unit cell lengths for a triclinic cell and will be optimised. Since there 
1420: are combinations of unit cell angles that are physically impossible, steps will be rejected if such combinations are detected. In \textit{mode} 2, as 
1421: well as the cell being optimised, potential parameters specified by \textit{params} (an integer between 1 and 3) 
1422: will also be optimised. If the 1 bit of \textit{params} is set, \textit{k} will be optimised, and if the 2 bit is set, $\phi$ will be optimised. 
1423:  
1424: \item {\it ORBITALS \/}: orbital localisation rotation potential.1406: \item {\it ORBITALS \/}: orbital localisation rotation potential.
1425: 1407: 
1426: \item {\it ORGYR \/}: radius of gyration will be calculated as an order parameter---requires 1408: \item {\it ORGYR \/}: radius of gyration will be calculated as an order parameter---requires 
1427: further documentation.1409: further documentation.
1428: 1410: 
1429: \item {\it OSASA \/}: order parameter---requires documentation.1411: \item {\it OSASA \/}: order parameter---requires documentation.
1430: 1412: 
1431: \item {\it P46\/}: specifies a 46-bead three-colour model polypeptide.1413: \item {\it P46\/}: specifies a 46-bead three-colour model polypeptide.
1432: See also the {\it BLN} keyword, which implements this potential in a more1414: See also the {\it BLN} keyword, which implements this potential in a more
1433: general way and uses unit bond lengths.1415: general way and uses unit bond lengths.
1841: requested for the run that produced the dumpfile, minus the number that were completed1823: requested for the run that produced the dumpfile, minus the number that were completed
1842: at the point the dumpfile was created. This option is not available before version 2.3.1824: at the point the dumpfile was created. This option is not available before version 2.3.
1843: If you are using the {\it A9INTE\/} keyword, you can specify the interaction enthalpy1825: If you are using the {\it A9INTE\/} keyword, you can specify the interaction enthalpy
1844: dump file to restore from as a second arguement.1826: dump file to restore from as a second arguement.
1845: 1827: 
1846: \item {\it RGCL2\/}: specifies a DIM rare gas-Cl$_2$ potential.1828: \item {\it RGCL2\/}: specifies a DIM rare gas-Cl$_2$ potential.
1847: 1829: 
1848: \item {\it RIGIDINIT\/}: This will turn on the local rigid body framework. This keyword needs two additional input files: {\it coordsinirigid\/} and {\it rbodyconfig\/}. {\it coordsinirigid\/} should contain the coordinates of all the atoms whether they are part of local rigid bodies or not. During initialization, GMIN will select appropriate coordinates and throw away irrelevant ones. The format is: \\1830: \item {\it RIGIDINIT\/}: This will turn on the local rigid body framework. This keyword needs two additional input files: {\it coordsinirigid\/} and {\it rbodyconfig\/}. {\it coordsinirigid\/} should contain the coordinates of all the atoms whether they are part of local rigid bodies or not. During initialization, GMIN will select appropriate coordinates and throw away irrelevant ones. The format is: \\
1849: {\it x1 y1 z1} \\1831: {\it x1 y1 z1} \\
1850: {\it x2 y2 z2} \\1832: {\it x2 y2 z2} \\
1851: {\it \ldots} \\1833: {\it ldots} \\
1852: {\it \ldots} \\1834: {\it ldots} \\
1853: {\it xn yn zn} \\1835: {\it xn yn zn} \\
1854: 1836: 
1855: The file {\it rbodyconfig\/} defines the rigid bodies, with the following format: \\1837: The file {\it rbodyconfig\/} defines the rigid bodies, with the following format: \\
1856: {\it GROUP NoAtomsInGroup} \\1838: {\it GROUP NoAtomsInGroup} \\
1857: {\it Atom 1} \\1839: {\it Atom 1} \\
1858: {\it Atom 2} \\1840: {\it Atom 2} \\
1859: {\it ldots} \\1841: {\it ldots} \\
1860: {\it ldots} \\1842: {\it ldots} \\
1861: {\it Atom N} \\1843: {\it Atom N} \\
1862: The keywords {\it RELAXFINALQUENCH\/} and {\it NRELAXRIGID\/} can be used in conjuction with {\it RIGIDINIT\/}.1844: The keywords {\it RELAXFINALQUENCH\/} and {\it NRELAXRIGID\/} can be used in conjuction with {\it RIGIDINIT\/}.
1960: displacements will be made for rigid bodies. Omitting {\it block\/} or using a value of zero results in1942: displacements will be made for rigid bodies. Omitting {\it block\/} or using a value of zero results in
1961: translational and orientational steps being taken simultaneously1943: translational and orientational steps being taken simultaneously
1962: for rigid bodies. The default values for {\it step\/},1944: for rigid bodies. The default values for {\it step\/},
1963: {\it astep\/} and {\it ostep\/} are all 0.3 and the default value of {\it nblock\/} is zero.1945: {\it astep\/} and {\it ostep\/} are all 0.3 and the default value of {\it nblock\/} is zero.
1964: 1946: 
1965: \item {\it STEEREDMIN smink sminkinc smindiststart smindistfinish sminatoma sminatomb\/}: specified steered 1947: \item {\it STEEREDMIN smink sminkinc smindiststart smindistfinish sminatoma sminatomb\/}: specified steered 
1966: minimisation should be performed (must be used with {\it AMBER9}). For a protein/ligand system, this adds a translation1948: minimisation should be performed (must be used with {\it AMBER9}). For a protein/ligand system, this adds a translation
1967: to the MC move. The vector between the centre of coordinates of groups A and B (as defined in the file movableatoms)1949: to the MC move. The vector between the centre of coordinates of groups A and B (as defined in the file movableatoms)
1968: is calculated and set to {\it smindiststart}. During the following minimisation, a restoring force is applied to 1950: is calculated and set to {\it smindiststart}. During the following minimisation, a restoring force is applied to 
1969: the ligand. The harmonic force constant is initially zero, and rises by {\it sminkinc} every LBFGS step up to a1951: the ligand. The harmonic force constant is initially zero, and rises by {\it sminkinc} every LBFGS step up to a
1970: maximum of {\it smink}. The force is applied until the A-\textgreater B distance is less than {\it smindistfinish}.  1952: maximum of {\it smink}. The force is applied until the A->B distance is less than {\it smindistfinish}.  
1971: 1953: 
1972: \item {\it STEPS mcsteps tfac\/}: determines the length of the1954: \item {\it STEPS mcsteps tfac\/}: determines the length of the
1973: basin-hopping run through the integer {\it mcsteps\/} and the annealing protocol through1955: basin-hopping run through the integer {\it mcsteps\/} and the annealing protocol through
1974: the real variable {\it tfac\/}. The temperature is multiplied by {\it tfac\/}1956: the real variable {\it tfac\/}. The temperature is multiplied by {\it tfac\/}
1975: after every step in each run. 1957: after every step in each run. 
1976: 1958: 
1977: \item {\it STICKY nrbsites, sigma\/}: specifies a `sticky patch' potential with {\it nrbsites}1959: \item {\it STICKY nrbsites, sigma\/}: specifies a `sticky patch' potential with {\it nrbsites}
1978: sites in the rigid body reference and a value of {\it sigma} for the $\sigma$ parameter.1960: sites in the rigid body reference and a value of {\it sigma} for the $\sigma$ parameter.
1979: 1961: 
1980: \item {\it STOCK mu lambda}: specifies a Stockmeyer potential with parameters1962: \item {\it STOCK mu lambda}: specifies a Stockmeyer potential with parameters


r32795/keywords.f 2017-06-14 20:30:11.311921006 +0100 r32794/keywords.f 2017-06-14 20:30:14.363960988 +0100
 36:       USE modamber 36:       USE modamber
 37:       USE PORFUNCS 37:       USE PORFUNCS
 38:       USE MYGA_PARAMS 38:       USE MYGA_PARAMS
 39:       USE BGUPMOD 39:       USE BGUPMOD
 40:       USE GLJYMOD 40:       USE GLJYMOD
 41:       USE CHIRO_MODULE, ONLY: CHIRO_SIGMA, CHIRO_MU, CHIRO_GAMMA, CHIRO_L 41:       USE CHIRO_MODULE, ONLY: CHIRO_SIGMA, CHIRO_MU, CHIRO_GAMMA, CHIRO_L
 42:       USE CONVEX_POLYHEDRA_MODULE, ONLY: INITIALISE_POLYHEDRA, K_COMPRESS, K_OVERLAP 42:       USE CONVEX_POLYHEDRA_MODULE, ONLY: INITIALISE_POLYHEDRA, K_COMPRESS, K_OVERLAP
 43:       USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS_MODE, LJ_GAUSS_RCUT, LJ_GAUSS_EPS, 43:       USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS_MODE, LJ_GAUSS_RCUT, LJ_GAUSS_EPS,
 44:      &                        LJ_GAUSS_R0, LJ_GAUSS_SIGMASQ, LJ_GAUSS_PARAMS, 44:      &                        LJ_GAUSS_R0, LJ_GAUSS_SIGMASQ, LJ_GAUSS_PARAMS,
 45:      &                        LJ_GAUSS_INITIALISE 45:      &                        LJ_GAUSS_INITIALISE
 46:       USE OPP_MOD, ONLY: OPP_MODE, OPP_RCUT, OPP_K, OPP_PHI, OPP_PARAMS, 
 47:      &                   OPP_INITIALISE 
 48:       USE LJ_VAR_RAD_MOD, ONLY: LJ_VAR_RAD_INITIALISE 
 49:       USE MBPOLMOD, ONLY: MBPOLINIT 46:       USE MBPOLMOD, ONLY: MBPOLINIT
 50:       USE SWMOD, ONLY: SWINIT, MWINIT 47:       USE SWMOD, ONLY: SWINIT, MWINIT
 51:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_GET_COORDS 48:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_GET_COORDS
 52: !     &                                 AMBER12_ATOMSS, AMBER12_SETUP, 49: !     &                                 AMBER12_ATOMSS, AMBER12_SETUP,
 53: !     &                                 AMBER12_RESIDUES, 50: !     &                                 AMBER12_RESIDUES,
 54: !     &                                 POPULATE_ATOM_DATA 51: !     &                                 POPULATE_ATOM_DATA
 55:       USE CHIRALITY, ONLY : CIS_TRANS_TOL 52:       USE CHIRALITY, ONLY : CIS_TRANS_TOL
 56:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR 53:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR
 57:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS, 54:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS,
 58:      &     PARSE_MLJ_PARAMS 55:      &     PARSE_MLJ_PARAMS
 97:       CHARACTER(LEN=100) TOPFILE,PARFILE 94:       CHARACTER(LEN=100) TOPFILE,PARFILE
 98:       CHARACTER(LEN=20) UNSTRING 95:       CHARACTER(LEN=20) UNSTRING
 99:       DOUBLE PRECISION LJREPBB, LJATTBB, LJREPLL, LJATTLL, LJREPNN, LJATTNN, 96:       DOUBLE PRECISION LJREPBB, LJATTBB, LJREPLL, LJATTLL, LJREPNN, LJATTNN,
100:      &                 HABLN, HBBLN, HCBLN, HDBLN, EABLN, EBBLN, ECBLN, EDBLN, TABLN, TBBLN, TCBLN, TDBLN 97:      &                 HABLN, HBBLN, HCBLN, HDBLN, EABLN, EBBLN, ECBLN, EDBLN, TABLN, TBBLN, TCBLN, TDBLN
101:       DOUBLE PRECISION LJREPBL, LJATTBL, LJREPBN, LJATTBN, LJREPLN, LJATTLN 98:       DOUBLE PRECISION LJREPBL, LJATTBL, LJREPBN, LJATTBN, LJREPLN, LJATTLN
102:  99: 
103: !     DC430 >100: !     DC430 >
104:       DOUBLE PRECISION :: LPL, LPR101:       DOUBLE PRECISION :: LPL, LPR
105:       LOGICAL          :: RBSYMTEST     ! jdf43>102:       LOGICAL          :: RBSYMTEST     ! jdf43>
106: 103: 
107: !      DOUBLE PRECISION :: VOL ! dj337104:       DOUBLE PRECISION :: VOL ! dj337
108: !105: !
109: !       sf344> added stuff106: !       sf344> added stuff
110: !107: !
111:       CHARACTER(LEN=10) check1108:       CHARACTER(LEN=10) check1
112:       CHARACTER(LEN=1) readswitch109:       CHARACTER(LEN=1) readswitch
113:       CHARACTER(LEN=4) J1CHAR110:       CHARACTER(LEN=4) J1CHAR
114:       CHARACTER(LEN=20) J2CHAR111:       CHARACTER(LEN=20) J2CHAR
115:       INTEGER iostatus, groupsize, groupatom,groupoffset,axis1,axis2,EOF112:       INTEGER iostatus, groupsize, groupatom,groupoffset,axis1,axis2,EOF
116:       INTEGER LUNIT, GETUNIT113:       INTEGER LUNIT, GETUNIT
117: 114: 
804: 801: 
805:       CAPSID=.FALSE.802:       CAPSID=.FALSE.
806:       STRANDT=.FALSE.803:       STRANDT=.FALSE.
807:       PAHT=.FALSE.804:       PAHT=.FALSE.
808:       TIP=.FALSE.805:       TIP=.FALSE.
809:       TTM3T=.FALSE.806:       TTM3T=.FALSE.
810:       QUADT=.FALSE.807:       QUADT=.FALSE.
811:       STOCKT=.FALSE.808:       STOCKT=.FALSE.
812:       LJCOULT=.FALSE.809:       LJCOULT=.FALSE.
813:       LJ_GAUSST=.FALSE.810:       LJ_GAUSST=.FALSE.
814:       OPPT=.FALSE. 
815:       COULN=0811:       COULN=0
816:       COULQ=0.0D0812:       COULQ=0.0D0
817:       COULSWAP = 0.0D0813:       COULSWAP = 0.0D0
818:       COULTEMP = 0.0D0814:       COULTEMP = 0.0D0
819:       GAYBERNET=.FALSE.815:       GAYBERNET=.FALSE.
820:       ELLIPSOIDT=.FALSE.816:       ELLIPSOIDT=.FALSE.
821:       PYGPERIODICT=.FALSE.817:       PYGPERIODICT=.FALSE.
822:       LJCAPSIDT=.FALSE.818:       LJCAPSIDT=.FALSE.
823:       PYBINARYT=.FALSE.819:       PYBINARYT=.FALSE.
824:       pyt = .false.820:       pyt = .false.
1125:       ORBITTOL=1.0D-31121:       ORBITTOL=1.0D-3
1126:       NOINVERSION=.FALSE.1122:       NOINVERSION=.FALSE.
1127: !1123: !
1128: ! ds656> Parallelised generalised basin-hopping1124: ! ds656> Parallelised generalised basin-hopping
1129:       GBHT=.FALSE.1125:       GBHT=.FALSE.
1130: 1126: 
1131: ! General mixed LJ systems1127: ! General mixed LJ systems
1132:       GLJT=.FALSE.1128:       GLJT=.FALSE.
1133: ! ds656> Multicomponent LJ system (different implementation to GLJ!)1129: ! ds656> Multicomponent LJ system (different implementation to GLJ!)
1134:       MLJT=.FALSE.1130:       MLJT=.FALSE.
1135: ! jwrm2> LJ with each atom having a different radius. Works with RIGIDINIT. 
1136:       LJ_VAR_RADT = .FALSE. 
1137: 1131: 
1138: ! khs26> Free energy basin-hopping stuff1132: ! khs26> Free energy basin-hopping stuff
1139:       FEBHT = .FALSE.1133:       FEBHT = .FALSE.
1140:       FETEMP = 0.0D01134:       FETEMP = 0.0D0
1141: ! khs26> This requires a minimum separation between the zero and non-zero1135: ! khs26> This requires a minimum separation between the zero and non-zero
1142: ! eigenvalues for runs using free energy basin-hopping. This is based on1136: ! eigenvalues for runs using free energy basin-hopping. This is based on
1143: ! the minimum found in quench zero.1137: ! the minimum found in quench zero.
1144:       MIN_ZERO_SEP = 0.0D01138:       MIN_ZERO_SEP = 0.0D0
1145:       MAX_ATTEMPTS = 51139:       MAX_ATTEMPTS = 5
1146: ! Use sparse hessian methods1140: ! Use sparse hessian methods
4787: 4781: 
4788: !4782: !
4789: ! NOT DOCUMENTED4783: ! NOT DOCUMENTED
4790: !4784: !
4791:       ELSE IF (WORD.EQ.'LJMF') THEN4785:       ELSE IF (WORD.EQ.'LJMF') THEN
4792:          LJMFT=.TRUE.4786:          LJMFT=.TRUE.
4793: !        CALL LJPARAMMF4787: !        CALL LJPARAMMF
4794:          WRITE(MYUNIT,'(A)') 'LJMF not currently maintained'4788:          WRITE(MYUNIT,'(A)') 'LJMF not currently maintained'
4795:          STOP4789:          STOP
4796: 4790: 
4797: ! jwrm2> LJ with each atom having a different radius. Works with RIGIDINIT. 
4798:       ELSE IF (WORD .EQ. 'LJVARRAD') THEN 
4799:          LJ_VAR_RADT = .TRUE. 
4800:          CALL LJ_VAR_RAD_INITIALISE() 
4801:  
4802: ! start = an N-oligomer is constructed by relocating NREPEAT units and placing them4791: ! start = an N-oligomer is constructed by relocating NREPEAT units and placing them
4803: ! at random distance R with Rmin <= R <= Rmax and angle alpha in the xy-plane4792: ! at random distance R with Rmin <= R <= Rmax and angle alpha in the xy-plane
4804: ! transxy: rigid body translation only in the xy-plane4793: ! transxy: rigid body translation only in the xy-plane
4805: ! rotz:  rigid body rotation only around the z-axis4794: ! rotz:  rigid body rotation only around the z-axis
4806: ! dihesc: only perturbation to the side chains4795: ! dihesc: only perturbation to the side chains
4807: !4796: !
4808:       ELSE IF (WORD.EQ.'MAKEOLIGO') THEN4797:       ELSE IF (WORD.EQ.'MAKEOLIGO') THEN
4809:          MAKEOLIGOT=.TRUE.4798:          MAKEOLIGOT=.TRUE.
4810:          CALL READA(UNSTRING)4799:          CALL READA(UNSTRING)
4811:          IF (TRIM(ADJUSTL(UNSTRING)).EQ.'START') THEN4800:          IF (TRIM(ADJUSTL(UNSTRING)).EQ.'START') THEN
6381: ! jwrm2>6370: ! jwrm2>
6382:       ELSE IF (WORD.EQ.'LJGAUSS') THEN6371:       ELSE IF (WORD.EQ.'LJGAUSS') THEN
6383:           LJ_GAUSST = .TRUE.6372:           LJ_GAUSST = .TRUE.
6384:           CALL READI(LJ_GAUSS_MODE)6373:           CALL READI(LJ_GAUSS_MODE)
6385:           CALL READF(LJ_GAUSS_RCUT)6374:           CALL READF(LJ_GAUSS_RCUT)
6386:           CALL READF(LJ_GAUSS_EPS)6375:           CALL READF(LJ_GAUSS_EPS)
6387:           CALL READF(LJ_GAUSS_R0)6376:           CALL READF(LJ_GAUSS_R0)
6388:           CALL READF(LJ_GAUSS_SIGMASQ)6377:           CALL READF(LJ_GAUSS_SIGMASQ)
6389:           IF (LJ_GAUSS_MODE .EQ. 3) CALL READI(LJ_GAUSS_PARAMS)6378:           IF (LJ_GAUSS_MODE .EQ. 3) CALL READI(LJ_GAUSS_PARAMS)
6390:           CALL LJ_GAUSS_INITIALISE()6379:           CALL LJ_GAUSS_INITIALISE()
6391:  
6392:       ELSE IF (WORD.EQ.'OPP') THEN 
6393:           OPPT = .TRUE. 
6394:           CALL READI(OPP_MODE) 
6395:           CALL READF(OPP_RCUT) 
6396:           CALL READF(OPP_K) 
6397:           CALL READF(OPP_PHI) 
6398:           IF (OPP_MODE .EQ. 2) CALL READI(OPP_PARAMS) 
6399:           CALL OPP_INITIALISE() 
6400: ! jwrm2> end6380: ! jwrm2> end
6401: 6381: 
6402:       ELSE IF (WORD.EQ.'STOCK') THEN6382:       ELSE IF (WORD.EQ.'STOCK') THEN
6403:          STOCKT=.TRUE.6383:          STOCKT=.TRUE.
6404:          RIGID=.TRUE.6384:          RIGID=.TRUE.
6405:          NRBSITES=16385:          NRBSITES=1
6406:          CALL READF(STOCKMU)6386:          CALL READF(STOCKMU)
6407:          CALL READF(STOCKLAMBDA)6387:          CALL READF(STOCKLAMBDA)
6408:          ALLOCATE(SITE(NRBSITES,3))6388:          ALLOCATE(SITE(NRBSITES,3))
6409: 6389: 


r32795/lj_gauss.f90 2017-06-14 20:30:12.343936294 +0100 r32794/lj_gauss.f90 2017-06-14 20:30:14.587963924 +0100
 61: !                  parameters LJ_GAUSS_EPS, LJ_GAUSS_R0, LJ_GAUSS_SIGMASQ and 61: !                  parameters LJ_GAUSS_EPS, LJ_GAUSS_R0, LJ_GAUSS_SIGMASQ and
 62: !                  the potential parameters will be optimised, depending on the 62: !                  the potential parameters will be optimised, depending on the
 63: !                  value of LJ_GAUSS_PARAMS 63: !                  value of LJ_GAUSS_PARAMS
 64: ! LJ_GAUSS_PARAMS: Specifies which of the three potential parameters will be 64: ! LJ_GAUSS_PARAMS: Specifies which of the three potential parameters will be
 65: !                  optimised. An integer between 1 and 7. 65: !                  optimised. An integer between 1 and 7.
 66: !                  1 bit set: LJ_GAUSS_EPS will be optimised 66: !                  1 bit set: LJ_GAUSS_EPS will be optimised
 67: !                  2 bit set: LJ_GAUSS_R0 will be optimised 67: !                  2 bit set: LJ_GAUSS_R0 will be optimised
 68: !                  3 bit set: LJ_GAUSS_SIGMASQ will be optimised 68: !                  3 bit set: LJ_GAUSS_SIGMASQ will be optimised
 69: INTEGER :: LJ_GAUSS_MODE, LJ_GAUSS_PARAMS 69: INTEGER :: LJ_GAUSS_MODE, LJ_GAUSS_PARAMS
 70:  70: 
 71: ! RCUTINV6, RCUTINV12: factors for LJ which don't depend on particle distance 71: ! RCUTINV6, RCUTINV12: factors for LJ which don't pened on particle distance
 72: DOUBLE PRECISION :: RCUTINV6, RCUTINV12 72: DOUBLE PRECISION :: RCUTINV6, RCUTINV12
 73:  73: 
 74: ! EXP_RCUT, PREF_RCUT: factors for Gauss which don't depend on the particle 74: ! EXP_RCUT, PREF_RCUT: factors for Gauss which don't depend on the particle
 75: !                      distance 75: !                      distance
 76: ! SCL_FCT: overall scale factor for the potential, to allow comparison between 76: ! SCL_FCT: overall scale factor for the potential, to allow comparison between
 77: !          different parameters 77: !          different parameters
 78: DOUBLE PRECISION :: EXP_RCUT, PREF_RCUT, SCL_FCT 78: DOUBLE PRECISION :: EXP_RCUT, PREF_RCUT, SCL_FCT
 79:  79: 
 80: ! PI: the usual meaning 80: ! PI: the usual meaning
 81: DOUBLE PRECISION, PARAMETER :: PI = 4*ATAN(1.D0) 81: DOUBLE PRECISION, PARAMETER :: PI = 4*ATAN(1.D0)
104: ! V_SIGMA: WCA cutoff distance104: ! V_SIGMA: WCA cutoff distance
105: ! x_LOWER, x_UPPER: bounds on potential parameters105: ! x_LOWER, x_UPPER: bounds on potential parameters
106: ! x_REPULSION: harmonic force constant for each potential parameter106: ! x_REPULSION: harmonic force constant for each potential parameter
107: DOUBLE PRECISION, PARAMETER :: V_EPS = 1.D-3107: DOUBLE PRECISION, PARAMETER :: V_EPS = 1.D-3
108: DOUBLE PRECISION, PARAMETER :: V_SIGMA = 3.D-1108: DOUBLE PRECISION, PARAMETER :: V_SIGMA = 3.D-1
109: DOUBLE PRECISION, PARAMETER :: EPS_LOWER = 0.D0, EPS_UPPER = 5.D0109: DOUBLE PRECISION, PARAMETER :: EPS_LOWER = 0.D0, EPS_UPPER = 5.D0
110: DOUBLE PRECISION, PARAMETER :: EPS_REPULSION = 1.D4110: DOUBLE PRECISION, PARAMETER :: EPS_REPULSION = 1.D4
111: DOUBLE PRECISION, PARAMETER :: R0_LOWER = 1.D0, R0_UPPER = 2.1D0111: DOUBLE PRECISION, PARAMETER :: R0_LOWER = 1.D0, R0_UPPER = 2.1D0
112: DOUBLE PRECISION, PARAMETER :: R0_REPULSION = 1.D4112: DOUBLE PRECISION, PARAMETER :: R0_REPULSION = 1.D4
113: DOUBLE PRECISION, PARAMETER :: SIGMASQ_LOWER = 0.01D0, SIGMASQ_UPPER = 0.03D0113: DOUBLE PRECISION, PARAMETER :: SIGMASQ_LOWER = 0.01D0, SIGMASQ_UPPER = 0.03D0
114: DOUBLE PRECISION, PARAMETER :: SIGMASQ_REPULSION = 1.D9114: DOUBLE PRECISION, PARAMETER :: SIGMASQ_REPULSION = 1.D6
115: 115: 
116: ! BOX_DERIV: stores information for calculating lattice derivatives116: ! BOX_DERIV: stores information for calculating lattice derivatives
117: TYPE :: BOX_DERIV117: TYPE :: BOX_DERIV
118:     ! ORTHOG: the orthogonal tansformation matrix, for changing from118:     ! ORTHOG: the orthogonal tansformation matrix, for changing from
119:     !         crystallographic space to orthogonal space first index is row,119:     !         crystallographic space to orthogonal space first index is row,
120:     !         second is column120:     !         second is column
121:     DOUBLE PRECISION :: ORTHOG(3,3)121:     DOUBLE PRECISION :: ORTHOG(3,3)
122: 122: 
123:     ! DERIV: derivatives of the orthogonal matrix wrt lattice parameters,123:     ! DERIV: derivatives of the orthogonal matrix wrt lattice parameters,
124:     !        first index indicates the row, second index the column and124:     !        first index indicates the row, second index the column and
472:         LJ12 = LJ6**2472:         LJ12 = LJ6**2
473:         ENERGY = LJ12 - 2.D0*LJ6473:         ENERGY = LJ12 - 2.D0*LJ6
474:         ! Correction to LJ to make potential go to zero at RCUT474:         ! Correction to LJ to make potential go to zero at RCUT
475:         ENERGY = ENERGY - (RCUTINV12 - 2*RCUTINV6)475:         ENERGY = ENERGY - (RCUTINV12 - 2*RCUTINV6)
476:         ! Correction to LJ to make derivative go to zero at RCUT476:         ! Correction to LJ to make derivative go to zero at RCUT
477:         ENERGY = ENERGY + (12.D0/LJ_GAUSS_RCUT) * &477:         ENERGY = ENERGY + (12.D0/LJ_GAUSS_RCUT) * &
478:                           (RCUTINV12 - RCUTINV6)* &478:                           (RCUTINV12 - RCUTINV6)* &
479:                           (MODRIJ - LJ_GAUSS_RCUT)479:                           (MODRIJ - LJ_GAUSS_RCUT)
480:         ! Correction to LJ to make second derivatives go to zero at480:         ! Correction to LJ to make second derivatives go to zero at
481:         ! RCUT481:         ! RCUT
482: !        ENERGY = ENERGY + (1.D0/LJ_GAUSS_RCUT**2) * &482:         ENERGY = ENERGY + (1.D0/LJ_GAUSS_RCUT**2) * &
483: !                          (42.D0*RCUTINV6 - 78.D0*RCUTINV12) * &483:                           (42.D0*RCUTINV6 - 78.D0*RCUTINV12) * &
484: !                          (MODRIJ - LJ_GAUSS_RCUT)**2484:                           (MODRIJ - LJ_GAUSS_RCUT)**2
485:         485:         
486:         ! Set output486:         ! Set output
487:         CALC_ENERGY_LJ = ENERGY487:         CALC_ENERGY_LJ = ENERGY
488: 488: 
489:     END FUNCTION CALC_ENERGY_LJ489:     END FUNCTION CALC_ENERGY_LJ
490: 490: 
491: !-------------------------------------------------------------------------------491: !-------------------------------------------------------------------------------
492: !492: !
493: ! Calculates the Gaussian contribution to the energy493: ! Calculates the Gaussian contribution to the energy
494: ! MODRIJ: distance beteen the particle pair494: ! MODRIJ: distance beteen the particle pair
512:         EXP_RIJ = LJ_GAUSS_EPS * &512:         EXP_RIJ = LJ_GAUSS_EPS * &
513:                   EXP(- (MODRIJ - LJ_GAUSS_R0)**2 / &513:                   EXP(- (MODRIJ - LJ_GAUSS_R0)**2 / &
514:                         (2.D0*LJ_GAUSS_SIGMASQ))514:                         (2.D0*LJ_GAUSS_SIGMASQ))
515:         ENERGY = - EXP_RIJ515:         ENERGY = - EXP_RIJ
516:         ! Correction to Gauss to make potential go to zero at RCUT516:         ! Correction to Gauss to make potential go to zero at RCUT
517:         ENERGY = ENERGY + EXP_RCUT517:         ENERGY = ENERGY + EXP_RCUT
518:         ! Correction to Gauss to make derivatives go to zero at RCUT518:         ! Correction to Gauss to make derivatives go to zero at RCUT
519:         ENERGY = ENERGY - PREF_RCUT*EXP_RCUT*(MODRIJ-LJ_GAUSS_RCUT)519:         ENERGY = ENERGY - PREF_RCUT*EXP_RCUT*(MODRIJ-LJ_GAUSS_RCUT)
520:         ! Correction to Gauss to make second derivatives go to zero at520:         ! Correction to Gauss to make second derivatives go to zero at
521:         ! RCUT521:         ! RCUT
522: !        ENERGY = ENERGY + EXP_RCUT * (PREF_RCUT**2 - 1.D0/LJ_GAUSS_SIGMASQ) * &522:         ENERGY = ENERGY + EXP_RCUT * (PREF_RCUT**2 - 1.D0/LJ_GAUSS_SIGMASQ) * &
523: !                          0.5D0*((MODRIJ - LJ_GAUSS_RCUT)**2)523:                           0.5D0*((MODRIJ - LJ_GAUSS_RCUT)**2)
524: 524: 
525:         ! Set output525:         ! Set output
526:         CALC_ENERGY_GAUSS = ENERGY526:         CALC_ENERGY_GAUSS = ENERGY
527: 527: 
528:     END FUNCTION CALC_ENERGY_GAUSS528:     END FUNCTION CALC_ENERGY_GAUSS
529: 529: 
530: !-------------------------------------------------------------------------------530: !-------------------------------------------------------------------------------
531: !531: !
532: ! Calculates the pair potential gradient for a pair of particles.532: ! Calculates the pair potential gradient for a pair of particles.
533: ! MODRIJ: distance beteen the particle pair533: ! MODRIJ: distance beteen the particle pair
553:         LJ6 = (1.D0/MODRIJ)**6553:         LJ6 = (1.D0/MODRIJ)**6
554:         LJ12 = LJ6**2554:         LJ12 = LJ6**2
555: 555: 
556:         ! Evaluate LJ derivatives556:         ! Evaluate LJ derivatives
557:         DVDR = DVDR + (12.D0/MODRIJ) * (LJ6 - LJ12)557:         DVDR = DVDR + (12.D0/MODRIJ) * (LJ6 - LJ12)
558:         ! Correction to LJ to make derivaitves go to zero at558:         ! Correction to LJ to make derivaitves go to zero at
559:         ! RCUT559:         ! RCUT
560:         DVDR = DVDR + (12.D0/LJ_GAUSS_RCUT) * (RCUTINV12 - RCUTINV6)560:         DVDR = DVDR + (12.D0/LJ_GAUSS_RCUT) * (RCUTINV12 - RCUTINV6)
561:         ! Correction to LJ to make second derivatives go to zero561:         ! Correction to LJ to make second derivatives go to zero
562:         ! at RCUT562:         ! at RCUT
563: !        DVDR = DVDR + (12.D0/LJ_GAUSS_RCUT**2) * &563:         DVDR = DVDR + (12.D0/LJ_GAUSS_RCUT**2) * &
564: !                      (7.D0*RCUTINV6 - 13.D0*RCUTINV12) * &564:                       (7.D0*RCUTINV6 - 13.D0*RCUTINV12) * &
565: !                      (MODRIJ - LJ_GAUSS_RCUT)565:                       (MODRIJ - LJ_GAUSS_RCUT)
566: 566: 
567:         ! Gauss factors567:         ! Gauss factors
568:         PREF_RIJ = (MODRIJ - LJ_GAUSS_R0) / LJ_GAUSS_SIGMASQ568:         PREF_RIJ = (MODRIJ - LJ_GAUSS_R0) / LJ_GAUSS_SIGMASQ
569:         EXP_RIJ = LJ_GAUSS_EPS * EXP(- (MODRIJ - LJ_GAUSS_R0)**2 / &569:         EXP_RIJ = LJ_GAUSS_EPS * EXP(- (MODRIJ - LJ_GAUSS_R0)**2 / &
570:                                        (2.D0*LJ_GAUSS_SIGMASQ))570:                                        (2.D0*LJ_GAUSS_SIGMASQ))
571: 571: 
572:         ! Evaluate Gauss derivatives572:         ! Evaluate Gauss derivatives
573:         DVDR = DVDR + PREF_RIJ * EXP_RIJ573:         DVDR = DVDR + PREF_RIJ * EXP_RIJ
574:         ! Correction to Gauss to make derivatives go to zero at574:         ! Correction to Gauss to make derivatives go to zero at
575:         ! RCUT575:         ! RCUT
576:         DVDR = DVDR - PREF_RCUT * EXP_RCUT576:         DVDR = DVDR - PREF_RCUT * EXP_RCUT
577:         ! Correction to Gauss to make second derivatives go to577:         ! Correction to Gauss to make second derivatives go to
578:         ! zero at RCUT578:         ! zero at RCUT
579: !        DVDR = DVDR - (EXP_RCUT / LJ_GAUSS_SIGMASQ - &579:         DVDR = DVDR - (EXP_RCUT / LJ_GAUSS_SIGMASQ - &
580: !                       (PREF_RCUT**2) * EXP_RCUT) * &580:                        (PREF_RCUT**2) * EXP_RCUT) * &
581: !                      (MODRIJ - LJ_GAUSS_RCUT)581:                       (MODRIJ - LJ_GAUSS_RCUT)
582: 582: 
583:         CALC_GRAD = DVDR583:         CALC_GRAD = DVDR
584: 584: 
585:     END FUNCTION CALC_GRAD585:     END FUNCTION CALC_GRAD
586: 586: 
587: !-------------------------------------------------------------------------------587: !-------------------------------------------------------------------------------
588: !588: !
589: ! Initialisation. Calculate some factors that won't be changing589: ! Initialisation. Calculate some factors that won't be changing
590: !590: !
591: !-------------------------------------------------------------------------------591: !-------------------------------------------------------------------------------
702:             END DO702:             END DO
703:         END IF703:         END IF
704: 704: 
705:         ! Box length steps705:         ! Box length steps
706:         IF (LJ_GAUSS_MODE .GE. 1) THEN706:         IF (LJ_GAUSS_MODE .GE. 1) THEN
707:             X(3*NATOMS-2:3*NATOMS) = X(3*NATOMS-2:3*NATOMS) + &707:             X(3*NATOMS-2:3*NATOMS) = X(3*NATOMS-2:3*NATOMS) + &
708:                                      VEC_RANDOM() * MAX_LENGTH_STEP708:                                      VEC_RANDOM() * MAX_LENGTH_STEP
709:         END IF709:         END IF
710: 710: 
711:         ! Box angle steps711:         ! Box angle steps
712:         IF (LJ_GAUSS_MODE .GE. 2) THEN712:         DO
713:             DO713:             IF (LJ_GAUSS_MODE .GE. 2) THEN
714:                 NEW_ANGLES(:) = X(3*NATOMS-5:3*NATOMS-3) + &714:                 NEW_ANGLES(:) = X(3*NATOMS-5:3*NATOMS-3) + &
715:                                 VEC_RANDOM() * MAX_ANGLE_STEP715:                                 VEC_RANDOM() * MAX_ANGLE_STEP
 716:             END IF
716: 717: 
717:             ! See whether the unit cell we've generated is valid718:             ! See whether the unit cell we've generated is valid
718:                 IF (CHECK_ANGLES(NEW_ANGLES(:))) THEN719:             IF (CHECK_ANGLES(NEW_ANGLES(:))) THEN
719:                     X(3*NATOMS-5:3*NATOMS-3) =  NEW_ANGLES(:)720:                 X(3*NATOMS-5:3*NATOMS-3) =  NEW_ANGLES(:)
720:                     EXIT721:                 EXIT
721:                 END IF722:             END IF
722:             END DO723:         END DO
723:         END IF 
724: 724: 
725:         ! Parameter steps725:         ! Parameter steps
726:         IF (LJ_GAUSS_MODE .EQ. 3) THEN726:         IF (LJ_GAUSS_MODE .EQ. 3) THEN
727:             IF (BTEST(LJ_GAUSS_PARAMS, 0)) THEN727:             IF (BTEST(LJ_GAUSS_PARAMS, 0)) THEN
728:                 CALL RANDOM_NUMBER(STEP_SIZE)728:                 CALL RANDOM_NUMBER(STEP_SIZE)
729:                 X(3*NATOMS-8) = X(3*NATOMS-8) + &729:                 X(3*NATOMS-8) = X(3*NATOMS-8) + &
730:                                 (STEP_SIZE - 0.5D0) * 2.D0 * MAX_EPS_STEP730:                                 (STEP_SIZE - 0.5D0) * 2.D0 * MAX_EPS_STEP
731:             END IF731:             END IF
732:             IF (BTEST(LJ_GAUSS_PARAMS, 1)) THEN732:             IF (BTEST(LJ_GAUSS_PARAMS, 1)) THEN
733:                 CALL RANDOM_NUMBER(STEP_SIZE)733:                 CALL RANDOM_NUMBER(STEP_SIZE)
1184:             ! because SCL_FCT does not depend on LJ_GAUSS_R0, but DVDR0 is more1184:             ! because SCL_FCT does not depend on LJ_GAUSS_R0, but DVDR0 is more
1185:             ! complicated than DVDE1185:             ! complicated than DVDE
1186:             ! Gaussian derivative1186:             ! Gaussian derivative
1187:             DVDR0 = - PREF_RIJ * EXP_RIJ1187:             DVDR0 = - PREF_RIJ * EXP_RIJ
1188:             ! Correction to make potential got to zero at RCUT1188:             ! Correction to make potential got to zero at RCUT
1189:             DVDR0 = DVDR0 + PREF_RCUT * EXP_RCUT1189:             DVDR0 = DVDR0 + PREF_RCUT * EXP_RCUT
1190:             ! Correction to make the derivate go to zero at RCUT1190:             ! Correction to make the derivate go to zero at RCUT
1191:             DVDR0 = DVDR0 + (1.D0 / LJ_GAUSS_SIGMASQ - PREF_RCUT**2) * &1191:             DVDR0 = DVDR0 + (1.D0 / LJ_GAUSS_SIGMASQ - PREF_RCUT**2) * &
1192:                              EXP_RCUT * (MODYIJ - LJ_GAUSS_RCUT)1192:                              EXP_RCUT * (MODYIJ - LJ_GAUSS_RCUT)
1193:             ! Correction to make the second derivative go to zero at RCUT1193:             ! Correction to make the second derivative go to zero at RCUT
1194: !            DVDR0 = DVDR0 - (3.D0 * PREF_RCUT/LJ_GAUSS_SIGMASQ - PREF_RCUT**3)*&1194:             DVDR0 = DVDR0 - (3.D0 * PREF_RCUT/LJ_GAUSS_SIGMASQ - PREF_RCUT**3)*&
1195: !                             EXP_RCUT * 0.5D0 * (MODYIJ - LJ_GAUSS_RCUT)**21195:                              EXP_RCUT * 0.5D0 * (MODYIJ - LJ_GAUSS_RCUT)**2
1196: 1196: 
1197:             ! Now we can do the LJ_GAUSS_R0 derivative1197:             ! Now we can do the LJ_GAUSS_R0 derivative
1198:             CALC_GRAD_PARAMS(2) = SCL_FCT * DVDR01198:             CALC_GRAD_PARAMS(2) = SCL_FCT * DVDR0
1199:         END IF1199:         END IF
1200: 1200: 
1201:         IF (BTEST(LJ_GAUSS_PARAMS, 2)) THEN1201:         IF (BTEST(LJ_GAUSS_PARAMS, 2)) THEN
1202:             ! Finally we do the derivative wrt LJ_GAUSS_SIGMASQ. This one is1202:             ! Finally we do the derivative wrt LJ_GAUSS_SIGMASQ. This one is
1203:             ! really complicated...1203:             ! really complicated...
1204:             ! Gaussian derivative1204:             ! Gaussian derivative
1205:             DVDS = - PREF_RIJ**2 * EXP_RIJ * 0.5D01205:             DVDS = - PREF_RIJ**2 * EXP_RIJ * 0.5D0
1206:             ! Correction to make potential got to zero at RCUT1206:             ! Correction to make potential got to zero at RCUT
1207:             DVDS = DVDS + PREF_RCUT**2 * EXP_RCUT * 0.5D01207:             DVDS = DVDS + PREF_RCUT**2 * EXP_RCUT * 0.5D0
1208:             ! Correction to make the derivative go to zero at RCUT1208:             ! Correction to make the derivative go to zero at RCUT
1209:             DVDS = DVDS + PREF_RCUT * EXP_RCUT * (MODYIJ - LJ_GAUSS_RCUT) * &1209:             DVDS = DVDS + PREF_RCUT * EXP_RCUT * (MODYIJ - LJ_GAUSS_RCUT) * &
1210:                    (1.D0 / LJ_GAUSS_SIGMASQ - PREF_RCUT**2 * 0.5D0)1210:                    (1.D0 / LJ_GAUSS_SIGMASQ - PREF_RCUT**2 * 0.5D0)
1211:             ! Correction to make the second derivative go to zero at RCUT1211:             ! Correction to make the second derivative go to zero at RCUT
1212: !            DVDS = DVDS + EXP_RCUT * 0.5D0 * (MODYIJ - LJ_GAUSS_RCUT)**2 * &1212:             DVDS = DVDS + EXP_RCUT * 0.5D0 * (MODYIJ - LJ_GAUSS_RCUT)**2 * &
1213: !                   (1.D0 / LJ_GAUSS_SIGMASQ**2 + PREF_RCUT**4 * 0.5D0 - &1213:                    (1.D0 / LJ_GAUSS_SIGMASQ**2 + PREF_RCUT**4 * 0.5D0 - &
1214: !                    5.D0 * PREF_RCUT**2 / (2.D0 * LJ_GAUSS_SIGMASQ))1214:                     5.D0 * PREF_RCUT**2 / (2.D0 * LJ_GAUSS_SIGMASQ))
1215: 1215: 
1216:             ! Derivative of SCL_FCT wrt LJ_GAUSS_SIGMASQ1216:             ! Derivative of SCL_FCT wrt LJ_GAUSS_SIGMASQ
1217:             DSFDS = - SCL_FCT**2 * PI * LJ_GAUSS_EPS / &1217:             DSFDS = - SCL_FCT**2 * PI * LJ_GAUSS_EPS / &
1218:                       SQRT(2.D0 * PI * LJ_GAUSS_SIGMASQ)1218:                       SQRT(2.D0 * PI * LJ_GAUSS_SIGMASQ)
1219: 1219: 
1220:             ! Now we have all the necessary factors1220:             ! Now we have all the necessary factors
1221:             CALC_GRAD_PARAMS(3) = SCL_FCT * DVDS + DSFDS * V1221:             CALC_GRAD_PARAMS(3) = SCL_FCT * DVDS + DSFDS * V
1222:         END IF1222:         END IF
1223: 1223: 
1224:     END FUNCTION CALC_GRAD_PARAMS1224:     END FUNCTION CALC_GRAD_PARAMS


r32795/lj_var_rad.f90 2017-06-14 20:30:12.563937407 +0100 r32794/lj_var_rad.f90 2017-06-14 20:30:14.807966807 +0100
  1: !   GMIN: A program for finding global minima  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/lj_var_rad.f90' in revision 32794
  2: !   Copyright (C) 1999-2017 David J. Wales 
  3: !   This file is part of GMIN. 
  4: ! 
  5: !   GMIN is free software; you can redistribute it and/or modify 
  6: !   it under the terms of the GNU General Public License as published by 
  7: !   the Free Software Foundation; either version 2 of the License, or 
  8: !   (at your option) any later version. 
  9: ! 
 10: !   GMIN is distributed in the hope that it will be useful, 
 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of 
 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 13: !   GNU General Public License for more details. 
 14: ! 
 15: !   You should have received a copy of the GNU General Public License 
 16: !   along with this program; if not, write to the Free Software 
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
 18: ! 
 19: ! 
 20: ! This module provides a routine for calculating the energy and gradients with 
 21: ! a double well Lennard-Jones + Gaussian potential. The possibility of varying 
 22: ! and optimising the potential parameneter is included. 
 23: ! Also included is a framework for a varying triclinic unit cell. This could be 
 24: ! separated out and generalised to all isotropic pair potentials. 
 25: ! Questions to jwrm2. 
 26: ! 
 27: MODULE LJ_VAR_RAD_MOD 
 28:   
 29: ! Public functions 
 30: PUBLIC :: LJ_VAR_RAD, LJ_VAR_RAD_INITIALISE 
 31: ! Private parameters 
 32: PRIVATE :: SPHERE_RADII, BODY_MAP, REJECT 
 33: ! Private functions 
 34: PRIVATE :: CALC_ENERGY, CALC_GRAD, CALC_BODY_MAP 
 35:  
 36: ! SPHERE_RADII: the radii of each of the LJ sphere, read in from a file 
 37: DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: SPHERE_RADII 
 38:  
 39: ! BODY_MAP: maps from the atomic index to the rigid body 
 40: INTEGER, ALLOCATABLE, DIMENSION(:) :: BODY_MAP 
 41:  
 42: CONTAINS 
 43:  
 44: !------------------------------------------------------------------------------- 
 45: ! 
 46: ! Main routine for calculating the energy and gradients. 
 47: ! X: coordinates array 
 48: ! GRAD: gradients array 
 49: ! ENERGY: calculated energy 
 50: ! GTEST: whether gradients should be calculated 
 51: ! 
 52: !------------------------------------------------------------------------------- 
 53:     SUBROUTINE LJ_VAR_RAD(X, GRAD, ENERGY, GTEST) 
 54:  
 55:         ! NATOMS: number of particles 
 56:         USE COMMONS, ONLY: NATOMS 
 57:          
 58:         IMPLICIT NONE 
 59:  
 60:         DOUBLE PRECISION, INTENT(INOUT) :: X(3*NATOMS) 
 61:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOMS), ENERGY 
 62:         LOGICAL, INTENT(IN) :: GTEST 
 63:  
 64:         ! I, J: indices for looping over particles 
 65:         ! NMOL: actual number of particles 
 66:         INTEGER :: I, J, NMOL 
 67:  
 68:         ! RIJ: Interparticle vector 
 69:         ! MODRIJ: distance between a pair of particles 
 70:         ! DVDR: derivative of pair potential wrt MODRIJ 
 71:         DOUBLE PRECISION :: RIJ(3), MODRIJ, DVDR 
 72:  
 73:         ! Allocate the BODY_MAP if it hasn't been done already. This will only 
 74:         ! be done once, but we can't do it with the other initialisation, as we 
 75:         ! can't be sure that the user has put RIGIDINIT before LJVARRAD in the 
 76:         ! data file. 
 77:         IF (.NOT. ALLOCATED(BODY_MAP)) CALL CALC_BODY_MAP(NATOMS) 
 78:  
 79:         ! Initialise output variables 
 80:         ENERGY = 0.D0 
 81:         GRAD(:) = 0.D0 
 82:  
 83:         NMOL = NATOMS 
 84:         DO I = 1, NMOL-1 ! Outer loop over particles 
 85:             DO J = I+1, NMOL ! Inner loop over particles 
 86:  
 87:                 ! Check whether atoms are in the same rigid body 
 88:                 IF (BODY_MAP(I) .NE. BODY_MAP(J)) THEN 
 89:                     ! Get the particle-particle distance 
 90:                     RIJ(:) = X(3*I-2:3*I) - X(3*J-2:3*J) 
 91:                     MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:))) 
 92:  
 93:                     ! There seems to be a problem with RIGIDINIT of some  
 94:                     ! particles ending up precisely on top of each other. Check 
 95:                     ! for this and reject. 
 96:                     IF (MODRIJ .LT. 1.D-03) THEN 
 97:                         CALL REJECT(ENERGY, GRAD(:), NATOMS) 
 98:                         RETURN 
 99:                     END IF 
100:  
101:                     ! Add energies and gradients 
102:                     ENERGY = ENERGY + CALC_ENERGY(MODRIJ, I, J) 
103:  
104:                     IF (GTEST) THEN 
105:                         DVDR = CALC_GRAD(MODRIJ, I, J) 
106:                         GRAD(3*I-2:3*I) = GRAD(3*I-2:3*I) + & 
107:                                           DVDR*RIJ(:)/MODRIJ 
108:                         GRAD(3*J-2:3*J) = GRAD(3*J-2:3*J) - & 
109:                                           DVDR*RIJ(:)/MODRIJ 
110: !                        WRITE (*, *) 'I = ', I, '; J = ', J, '; DVDR = ', DVDR, '; MODRIJ = ', MODRIJ 
111:                     END IF ! GTEST 
112:                 END IF ! If less than cutoff distance 
113:             END DO ! Inner loop over particles 
114:         END DO ! Outer loop over particles 
115:  
116: !    IF (GTEST) WRITE (*, *) 'LJ_VAR_RAD> GRAD = ', GRAD 
117: !    WRITE (MYUNIT, *) 'LJ_VAR_RAD> NOMR2(GRAD) = ', NORM2(GRAD) 
118:  
119:     END SUBROUTINE LJ_VAR_RAD 
120:  
121: !------------------------------------------------------------------------------- 
122: ! 
123: ! Rejects a step by setting the energy very high and the gradients very low. 
124: ! This tricks L-BFGS into thinking it's converged, but the high energy of the 
125: ! quench will ensure it is rejected. 
126: ! ENERGY: system energy 
127: ! GRAD: system gradients 
128: ! NATOMS: number of atoms in the system 
129: ! 
130: !------------------------------------------------------------------------------- 
131:  
132:     SUBROUTINE REJECT(ENERGY, GRAD, NATOMS) 
133:  
134:         IMPLICIT NONE 
135:  
136:         INTEGER, INTENT(IN) :: NATOMS 
137:         DOUBLE PRECISION, INTENT(OUT) :: ENERGY, GRAD(3*NATOMS) 
138:  
139:         ENERGY = 1.D20 
140:         GRAD(:) = 1.D-20 
141:  
142:     END SUBROUTINE REJECT 
143:  
144: !------------------------------------------------------------------------------- 
145: ! 
146: ! Calculates the LJ contribution to the energy 
147: ! MODRIJ: distance beteen the particle pair 
148: ! I, J: indices of atoms 
149: ! 
150: !------------------------------------------------------------------------------- 
151:  
152:     PURE DOUBLE PRECISION FUNCTION CALC_ENERGY (MODRIJ, I, J) 
153:  
154:         IMPLICIT NONE 
155:  
156:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
157:         INTEGER, INTENT(IN) :: I, J 
158:  
159:         ! DELTA: adjustment to radius for this pair 
160:         DOUBLE PRECISION :: DELTA 
161:  
162:         DELTA = 0.5D0 * (SPHERE_RADII(I) + SPHERE_RADII(J)) - 1.D0 
163:         CALC_ENERGY = 4.D0 * ((1.D0 / (MODRIJ - DELTA))**12 - & 
164:                               (1.D0 / (MODRIJ - DELTA))**6) 
165:  
166:     END FUNCTION CALC_ENERGY 
167:  
168: !------------------------------------------------------------------------------- 
169: ! 
170: ! Calculates the pair potential gradient for a pair of particles. 
171: ! MODRIJ: distance beteen the particle pair 
172: ! I, J: indices of atoms 
173: ! 
174: !------------------------------------------------------------------------------- 
175:     PURE DOUBLE PRECISION FUNCTION CALC_GRAD (MODRIJ, I, J) 
176:  
177:         IMPLICIT NONE 
178:  
179:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
180:         INTEGER, INTENT(IN) :: I, J 
181:  
182:         ! DELTA: adjustment to radius for this pair 
183:         DOUBLE PRECISION :: DELTA 
184:  
185:         DELTA = 0.5D0 * (SPHERE_RADII(I) + SPHERE_RADII(J)) - 1.D0 
186:         CALC_GRAD = 24.D0 * ((1.D0 / (MODRIJ - DELTA))**7 - & 
187:                       2.D0 * (1.D0 / (MODRIJ - DELTA))**13) 
188:  
189:     END FUNCTION CALC_GRAD 
190:  
191: !------------------------------------------------------------------------------- 
192: ! 
193: ! Initialisation. Read in the radii from the file 'lj_var_rad' 
194: ! 
195: !------------------------------------------------------------------------------- 
196:     SUBROUTINE LJ_VAR_RAD_INITIALISE () 
197:  
198:         ! MYUNIT: file unit for main GMIN output 
199:         ! NATOMS: number of atoms in the system 
200:         USE COMMONS, ONLY: MYUNIT, NATOMS 
201:  
202:         IMPLICIT NONE 
203:  
204:         ! GETUNIT: function for fiding a free file unit 
205:         ! RADIUS_UNIT: file unit for the radii file 
206:         ! OPEN_STATUS: for checking opening of the file 
207:         ! I: loop variable 
208:         INTEGER :: GETUNIT, RADIUS_UNIT, OPEN_STATUS, I 
209:  
210:         ! Open the file, checking it already exists 
211:         RADIUS_UNIT = GETUNIT() 
212:         OPEN(UNIT=RADIUS_UNIT, FILE='lj_var_rad', STATUS='OLD', &    
213:              IOSTAT=OPEN_STATUS, ACTION='READ') 
214:         IF (OPEN_STATUS .NE. 0) THEN 
215:             WRITE(MYUNIT, *) & 
216:                 'LJ_VAR_RAD_INITIALISE> Could not open lj_var_rad for reading' 
217:             STOP 
218:         END IF 
219:  
220:         ALLOCATE(SPHERE_RADII(NATOMS)) 
221:         SPHERE_RADII(:) = 1.D0 
222:         DO I = 1, NATOMS 
223:             READ(RADIUS_UNIT, *) SPHERE_RADII(I) 
224:         END DO 
225:  
226:         CLOSE(RADIUS_UNIT) 
227:  
228:     END SUBROUTINE LJ_VAR_RAD_INITIALISE 
229:  
230: !------------------------------------------------------------------------------- 
231: ! 
232: ! Initialisation. Work out the mapping from atomic index to rigid body index 
233: ! NATOMS: number of atoms in the system 
234: ! 
235: !------------------------------------------------------------------------------- 
236:     SUBROUTINE CALC_BODY_MAP(NATOMS) 
237:  
238:         ! NRIGIDBODIES: the number of rigid bodies 
239:         ! NSITEPERBODY: the number of atoms for each body 
240:         ! RIGIDINIT: whether the rigid body framework is active 
241:         USE GENRIGID, ONLY: NRIGIDBODY, NSITEPERBODY, RIGIDINIT 
242:  
243:         IMPLICIT NONE 
244:  
245:         INTEGER, INTENT(IN) :: NATOMS 
246:  
247:         ! I, J: loop indices 
248:         INTEGER :: I, J, K 
249:  
250:         ALLOCATE(BODY_MAP(NATOMS)) 
251:  
252:         IF (RIGIDINIT) THEN 
253:             ! K will be our atom counter 
254:             K = 1 
255:             DO I = 1, NRIGIDBODY ! loop over each body 
256:                 DO J = 1, NSITEPERBODY(I) ! loop over the atoms within the body 
257:                     BODY_MAP(K) = I 
258:                     K = K + 1 
259:                 END DO 
260:             END DO 
261:         ELSE ! Not using rigid body framework 
262:             DO I = 1, NATOMS 
263:                 BODY_MAP(I) = I 
264:             END DO 
265:         END IF 
266:  
267:     END SUBROUTINE CALC_BODY_MAP 
268:  
269: END MODULE LJ_VAR_RAD_MOD 


r32795/mc.F 2017-06-14 20:30:12.787940340 +0100 r32794/mc.F 2017-06-14 20:30:15.031969738 +0100
 30:      &                      SETCHIRAL,AMCHPMAX,DOLIGMOVE,LIGMOVEFREQ, AMCHNMAX, ROTAMERT, 30:      &                      SETCHIRAL,AMCHPMAX,DOLIGMOVE,LIGMOVEFREQ, AMCHNMAX, ROTAMERT,
 31:      &                      E_IGB, E_BOND, E_ANGLE, E_DIHEDRAL, 31:      &                      E_IGB, E_BOND, E_ANGLE, E_DIHEDRAL,
 32:      &                      E_VDW, E_14_VDW, E_ELEC, E_14_ELEC, IGB, MACROIONT 32:      &                      E_VDW, E_14_VDW, E_ELEC, E_14_ELEC, IGB, MACROIONT
 33:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_WRITE_RESTART, AMBER12_WRITE_PDB, POT_ENE_REC_C 33:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_WRITE_RESTART, AMBER12_WRITE_PDB, POT_ENE_REC_C
 34:       USE CHIRALITY, ONLY : INIT_CHIRAL, INIT_CIS_TRANS 34:       USE CHIRALITY, ONLY : INIT_CHIRAL, INIT_CIS_TRANS
 35:       USE porfuncs 35:       USE porfuncs
 36:       USE AMHGLOBALS, ONLY: NMRES,OMOVI,AVEP,NUMPRO,IRES 36:       USE AMHGLOBALS, ONLY: NMRES,OMOVI,AVEP,NUMPRO,IRES
 37:       USE AMH_INTERFACES, ONLY:E_WRITE 37:       USE AMH_INTERFACES, ONLY:E_WRITE
 38:       USE ROTAMER 38:       USE ROTAMER
 39:       USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS_TAKESTEP 39:       USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS_TAKESTEP
 40:       USE OPP_MOD, ONLY: OPP_TAKESTEP 
 41:  40: 
 42:       IMPLICIT NONE 41:       IMPLICIT NONE
 43: #ifdef MPI 42: #ifdef MPI
 44:       INCLUDE 'mpif.h' 43:       INCLUDE 'mpif.h'
 45:       INTEGER MPIERR 44:       INTEGER MPIERR
 46:       LOGICAL HITANY,MPI_FINT 45:       LOGICAL HITANY,MPI_FINT
 47: #endif 46: #endif
 48:  47: 
 49:       INTEGER J1, NSUCCESS(NPAR), NFAIL(NPAR), NFAILT(NPAR), NSUCCESST(NPAR), J2, NSTEPS, JP, REJSTREAK, 48:       INTEGER J1, NSUCCESS(NPAR), NFAIL(NPAR), NFAILT(NPAR), NSUCCESST(NPAR), J2, NSTEPS, JP, REJSTREAK,
 50:      1        UNT, ITERATIONS, NSUPERCOUNT, NQTOT, JACCPREV, NREN, NLAST, NSTEPREN, BRUN,QDONE,JBEST(NPAR),GCJBEST(NPAR), 49:      1        UNT, ITERATIONS, NSUPERCOUNT, NQTOT, JACCPREV, NREN, NLAST, NSTEPREN, BRUN,QDONE,JBEST(NPAR),GCJBEST(NPAR),
 84:       COMMON /Q4C/ QSTART, QFINISH 83:       COMMON /Q4C/ QSTART, QFINISH
 85:  84: 
 86:       character(len=10)       :: datechar,timechar,zonechar 85:       character(len=10)       :: datechar,timechar,zonechar
 87:       integer                 :: values(8),itime1 86:       integer                 :: values(8),itime1
 88:       DOUBLE PRECISION :: DISTGROUPX2,DISTGROUPY2,DISTGROUPZ2,DISTGROUPCENTRE,DUMMY 87:       DOUBLE PRECISION :: DISTGROUPX2,DISTGROUPY2,DISTGROUPZ2,DISTGROUPCENTRE,DUMMY
 89:       integer :: J6, L, M, SUMSQUAREDIFF, SUMSQUAREDIFF1D, HBONDDIFF 88:       integer :: J6, L, M, SUMSQUAREDIFF, SUMSQUAREDIFF1D, HBONDDIFF
 90:       INTEGER GETUNIT 89:       INTEGER GETUNIT
 91:       INTEGER DUMPUNIQUEUNIT 90:       INTEGER DUMPUNIQUEUNIT
 92:  91: 
 93: ! khs26> Energy decomposition for AMBER 12 92: ! khs26> Energy decomposition for AMBER 12
 94: !      TYPE(POT_ENE_REC_C) :: AMBER12_ENERGY_DECOMP 93:       TYPE(POT_ENE_REC_C) :: AMBER12_ENERGY_DECOMP
 95:  94: 
 96: ! ds656 > energy before and after BGUPTA swap (without relaxation) 95: ! ds656 > energy before and after BGUPTA swap (without relaxation)
 97:       DOUBLE PRECISION :: EBSWAP=0.0D0, EASWAP=0.0D0, EASWAPQ=0.0D0, DE1=0.0D0, DE2=0.0D0 96:       DOUBLE PRECISION :: EBSWAP=0.0D0, EASWAP=0.0D0, EASWAPQ=0.0D0, DE1=0.0D0, DE2=0.0D0
 98:  97: 
 99:       DOUBLE PRECISION QTEMP(3,NATOMSALLOC), Q4, Q6 98:       DOUBLE PRECISION QTEMP(3,NATOMSALLOC), Q4, Q6
100:  99: 
101:       QNEWRES=0100:       QNEWRES=0
102:       QGCBH=0101:       QGCBH=0
103: 102: 
104: ! sf344> selectively unfreeze the last added building blocks (MODULARNINCR)103: ! sf344> selectively unfreeze the last added building blocks (MODULARNINCR)
985: 984: 
986: !-------------------------------------------985: !-------------------------------------------
987: ! hk286 - generalised THOMSON986: ! hk286 - generalised THOMSON
988: !-----------------------------------------------987: !-----------------------------------------------
989:                ELSE IF (GTHOMSONT) THEN988:                ELSE IF (GTHOMSONT) THEN
990:                 CALL TAKESTEPGTHOMSON ()989:                 CALL TAKESTEPGTHOMSON ()
991: 990: 
992:                ELSE IF (GAYBERNEDCT) THEN991:                ELSE IF (GAYBERNEDCT) THEN
993:                 CALL TAKESTEPGB(JP)992:                 CALL TAKESTEPGB(JP)
994: 993: 
995: ! jwrm2> We need a special step taking routine for the LJ_GAUSS and OPP, to deal994: ! jwrm2> We need a special step taking routine for the LJGAUSS, to deal
996: !        getting non-physical unit cell and potential parameters995: !        getting non-physical unit cell and potential parameters
997:                ELSE IF (LJ_GAUSST) THEN996:                ELSE IF (LJ_GAUSST) THEN
998:                 CALL LJ_GAUSS_TAKESTEP(JP)997:                 CALL LJ_GAUSS_TAKESTEP(JP)
999: 998: 
1000:                ELSE IF (OPPT) THEN 
1001:                 CALL OPP_TAKESTEP(JP) 
1002:  
1003:                ELSE IF (PYT) THEN999:                ELSE IF (PYT) THEN
1004: !               seed the random number generator with system time + MYNODE (for MPI runs)1000: !               seed the random number generator with system time + MYNODE (for MPI runs)
1005:                   IF(RANDOMSEEDT) THEN1001:                   IF(RANDOMSEEDT) THEN
1006:                      CALL DATE_AND_TIME(datechar,timechar,zonechar,values)1002:                      CALL DATE_AND_TIME(datechar,timechar,zonechar,values)
1007:                      itime1= values(6)*60 + values(7)1003:                      itime1= values(6)*60 + values(7)
1008:                      CALL SDPRND(itime1+MYNODE)1004:                      CALL SDPRND(itime1+MYNODE)
1009:                   END IF1005:                   END IF
1010:                 call py_takestep(jp)1006:                 call py_takestep(jp)
1011:                 if (rigidmdt) call rigidmd_nvt1007:                 if (rigidmdt) call rigidmd_nvt
1012: 1008: 
1230:                         SYMFCTR=SYMFCTR+LOG(1.0D0*NATOMS)1226:                         SYMFCTR=SYMFCTR+LOG(1.0D0*NATOMS)
1231:                         IF (MOD(J1-1,PRTFRQ).EQ.0) WRITE(MYUNIT,'(A,G20.10)') 1227:                         IF (MOD(J1-1,PRTFRQ).EQ.0) WRITE(MYUNIT,'(A,G20.10)') 
1232:      &  'mc> Attempting to add an extra atom at radial distance from centre of coordinates ',CMMAX1228:      &  'mc> Attempting to add an extra atom at radial distance from centre of coordinates ',CMMAX
1233: !1229: !
1234: ! To generate uniformly distributed random points on the (n − 1)-sphere (i.e., 1230: ! To generate uniformly distributed random points on the (n − 1)-sphere (i.e., 
1235: ! the surface of the n-ball), Marsaglia (1972) gives the following algorithm.1231: ! the surface of the n-ball), Marsaglia (1972) gives the following algorithm.
1236: ! Generate an n-dimensional vector of normal deviates (it suffices to use N(0, 1), although in fact the choice of the variance is arbitrary),1232: ! Generate an n-dimensional vector of normal deviates (it suffices to use N(0, 1), although in fact the choice of the variance is arbitrary),
1237: ! Now calculate the "radius" of this point. 1233: ! Now calculate the "radius" of this point. 
1238: ! The vector is uniformly distributed over the surface of the unit n-ball.1234: ! The vector is uniformly distributed over the surface of the unit n-ball.
1239: !1235: !
1240: ! Marsaglia polar method is a modification of the Box-Muller method algorithm, 1236: ! Marsaglia polar method is a modification of the BoxMuller method algorithm, 
1241: ! which does not require computation of functions sin() and cos(). 1237: ! which does not require computation of functions sin() and cos(). 
1242: ! In this method U and V are drawn from the uniform (-1,1) distribution, 1238: ! In this method U and V are drawn from the uniform (1,1) distribution, 
1243: ! and then S = U2 + V2 is computed. If S is greater or equal to one then the method starts over, otherwise two quantities1239: ! and then S = U2 + V2 is computed. If S is greater or equal to one then the method starts over, otherwise two quantities
1244: ! are returned. Again, X and Y will be independent and standard normally distributed.1240: ! are returned. Again, X and Y will be independent and standard normally distributed.
1245: !1241: !
1246: ! The Box-Muller method uses two independent random numbers U and V distributed uniformly on (0,1). Then the two random variables X and Y1242: ! The BoxMuller method uses two independent random numbers U and V distributed uniformly on (0,1). Then the two random variables X and Y
1247: ! will both have the standard normal distribution, and will be independent. This formulation arises because for a bivariate normal random 1243: ! will both have the standard normal distribution, and will be independent. This formulation arises because for a bivariate normal random 
1248: ! vector (X Y) the squared norm X^2 + Y^2 will have the chi-squared distribution with two degrees of freedom, which is an easily generated 1244: ! vector (X Y) the squared norm X^2 + Y^2 will have the chi-squared distribution with two degrees of freedom, which is an easily generated 
1249: ! exponential random variable corresponding to the quantity -2ln(U) in these equations; and the angle is distributed uniformly around the circle, 1245: ! exponential random variable corresponding to the quantity 2ln(U) in these equations; and the angle is distributed uniformly around the circle, 
1250: ! chosen by the random variable V.1246: ! chosen by the random variable V.
1251: !1247: !
1252: 567                     DUMMY1=DPRAND()1248: 567                     DUMMY1=DPRAND()
1253:                         DUMMY2=DPRAND()1249:                         DUMMY2=DPRAND()
1254:                         COORDS(3*NATOMS-2,JP)=SQRT(-2.0D0*LOG(DUMMY1))*COS(6.2831853071796D0*DUMMY2)1250:                         COORDS(3*NATOMS-2,JP)=SQRT(-2.0D0*LOG(DUMMY1))*COS(6.2831853071796D0*DUMMY2)
1255:                         COORDS(3*NATOMS-1,JP)=SQRT(-2.0D0*LOG(DUMMY1))*SIN(6.2831853071796D0*DUMMY2)1251:                         COORDS(3*NATOMS-1,JP)=SQRT(-2.0D0*LOG(DUMMY1))*SIN(6.2831853071796D0*DUMMY2)
1256:                         DUMMY1=DPRAND()1252:                         DUMMY1=DPRAND()
1257:                         DUMMY2=DPRAND()1253:                         DUMMY2=DPRAND()
1258:                         COORDS(3*NATOMS,JP)=SQRT(-2.0D0*LOG(DUMMY1))*COS(6.2831853071796D0*DUMMY2)1254:                         COORDS(3*NATOMS,JP)=SQRT(-2.0D0*LOG(DUMMY1))*COS(6.2831853071796D0*DUMMY2)
1259:                         DUMMY1=SQRT(COORDS(3*NATOMS-2,JP)**2+COORDS(3*NATOMS-1,JP)**2+COORDS(3*NATOMS,JP)**2)1255:                         DUMMY1=SQRT(COORDS(3*NATOMS-2,JP)**2+COORDS(3*NATOMS-1,JP)**2+COORDS(3*NATOMS,JP)**2)


r32795/modcudalbfgs.F90 2017-06-14 20:30:13.007943221 +0100 r32794/modcudalbfgs.F90 2017-06-14 20:30:15.255972675 +0100
281:     END SUBROUTINE CUDA_LBFGS_WRAPPER281:     END SUBROUTINE CUDA_LBFGS_WRAPPER
282: 282: 
283:     SUBROUTINE CUDA_ENEGRAD_WRAPPER(NATOMS, COORDS, TOTENERGY, GRADIENTS)283:     SUBROUTINE CUDA_ENEGRAD_WRAPPER(NATOMS, COORDS, TOTENERGY, GRADIENTS)
284: 284: 
285:         INTEGER(KIND=C_INT) :: NATOMS, C_NADDTARGET285:         INTEGER(KIND=C_INT) :: NATOMS, C_NADDTARGET
286: 286: 
287:         REAL(KIND=C_DOUBLE) :: C_TOTENERGY, POTENTIALTIME287:         REAL(KIND=C_DOUBLE) :: C_TOTENERGY, POTENTIALTIME
288: #ifndef DUMMY_CUDA288: #ifndef DUMMY_CUDA
289:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS) :: COORDS, C_GRADIENTS289:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS) :: COORDS, C_GRADIENTS
290:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT290:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT
291: #else /* ifdef DUMMY_CUDA */291: #else ifdef DUMMY_CUDA
292:         REAL(KIND=C_DOUBLE), DIMENSION(1) :: COORDS, C_GRADIENTS292:         REAL(KIND=C_DOUBLE), DIMENSION(1) :: COORDS, C_GRADIENTS
293:         REAL(KIND=C_DOUBLE), DIMENSION(1, 1) :: C_LJADDREP, C_LJADDATT293:         REAL(KIND=C_DOUBLE), DIMENSION(1, 1) :: C_LJADDREP, C_LJADDATT
294: #endif294: #endif
295: 295: 
296:         LOGICAL(KIND=C_BOOL) :: C_CUDATIMET296:         LOGICAL(KIND=C_BOOL) :: C_CUDATIMET
297: 297: 
298:         CHARACTER(LEN=1, KIND=C_CHAR) :: C_CUDAPOT298:         CHARACTER(LEN=1, KIND=C_CHAR) :: C_CUDAPOT
299: 299: 
300:         INTEGER :: X, I, J300:         INTEGER :: X, I, J
301: 301: 
333:     END SUBROUTINE CUDA_ENEGRAD_WRAPPER333:     END SUBROUTINE CUDA_ENEGRAD_WRAPPER
334: 334: 
335:     SUBROUTINE CUDA_NUMERICAL_HESS(NATOMS, COORDS, HESSIAN, DELTA)335:     SUBROUTINE CUDA_NUMERICAL_HESS(NATOMS, COORDS, HESSIAN, DELTA)
336:         IMPLICIT NONE336:         IMPLICIT NONE
337: 337: 
338:         INTEGER(KIND=C_INT) :: NATOMS, C_NADDTARGET338:         INTEGER(KIND=C_INT) :: NATOMS, C_NADDTARGET
339:         REAL(KIND=C_DOUBLE) :: C_ENERGY, POTENTIALTIME339:         REAL(KIND=C_DOUBLE) :: C_ENERGY, POTENTIALTIME
340:         REAL(KIND=C_DOUBLE) :: COORDS(3*NATOMS), C_GRADIENTS(3*NATOMS)340:         REAL(KIND=C_DOUBLE) :: COORDS(3*NATOMS), C_GRADIENTS(3*NATOMS)
341: #ifndef DUMMY_CUDA341: #ifndef DUMMY_CUDA
342:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT342:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT
343: #else /* ifdef DUMMY_CUDA */343: #else ifdef DUMMY_CUDA
344:         REAL(KIND=C_DOUBLE), DIMENSION(1, 1) :: C_LJADDREP, C_LJADDATT344:         REAL(KIND=C_DOUBLE), DIMENSION(1, 1) :: C_LJADDREP, C_LJADDATT
345: #endif345: #endif
346: 346: 
347:         LOGICAL(KIND=C_BOOL) :: C_CUDATIMET347:         LOGICAL(KIND=C_BOOL) :: C_CUDATIMET
348: 348: 
349:         CHARACTER(LEN=1, KIND=C_CHAR) :: C_CUDAPOT349:         CHARACTER(LEN=1, KIND=C_CHAR) :: C_CUDAPOT
350: 350: 
351:         DOUBLE PRECISION    :: HESSIAN(3*NATOMS, 3*NATOMS)351:         DOUBLE PRECISION    :: HESSIAN(3*NATOMS, 3*NATOMS)
352:         DOUBLE PRECISION    :: DELTA352:         DOUBLE PRECISION    :: DELTA
353:         DOUBLE PRECISION    :: GRAD_PLUS(3*NATOMS), GRAD_MINUS(3*NATOMS)353:         DOUBLE PRECISION    :: GRAD_PLUS(3*NATOMS), GRAD_MINUS(3*NATOMS)


r32795/opp.f90 2017-06-14 20:30:13.231946156 +0100 r32794/opp.f90 2017-06-14 20:30:15.471975506 +0100
  1: !   GMIN: A program for finding global minima  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/opp.f90' in revision 32794
  2: !   Copyright (C) 1999-2017 David J. Wales 
  3: !   This file is part of GMIN. 
  4: ! 
  5: !   GMIN is free software; you can redistribute it and/or modify 
  6: !   it under the terms of the GNU General Public License as published by 
  7: !   the Free Software Foundation; either version 2 of the License, or 
  8: !   (at your option) any later version. 
  9: ! 
 10: !   GMIN is distributed in the hope that it will be useful, 
 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of 
 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 13: !   GNU General Public License for more details. 
 14: ! 
 15: !   You should have received a copy of the GNU General Public License 
 16: !   along with this program; if not, write to the Free Software 
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
 18: ! 
 19: ! 
 20: ! This module provides a routine for calculating the energy and gradients with 
 21: ! an oscillating pair potential (OPP). The possibility of varying 
 22: ! and optimising the potential parameters is included. 
 23: ! Also included is a framework for a varying triclinic unit cell. This could be 
 24: ! separated out and generalised to all isotropic pair potentials. 
 25: ! Questions to jwrm2. 
 26: ! 
 27: MODULE OPP_MOD 
 28:  
 29: ! Public parameters 
 30: PUBLIC :: OPP_MODE, OPP_RCUT, OPP_K, OPP_PHI, OPP_PARAMS 
 31: ! Public functions 
 32: PUBLIC :: OPP, OPP_INITIALISE, VIEW_OPP, OPP_TAKESTEP 
 33: ! Private parameters 
 34: PRIVATE :: V_EPS, V_SIGMA, FD_RCUT, SD_RCUT, IMAGE_CUTOFF 
 35: PRIVATE :: MAX_LENGTH_STEP, MAX_ANGLE_STEP, MAX_K_STEP, MAX_PHI_STEP 
 36: PRIVATE :: K_REPULSION, K_LOWER, K_UPPER 
 37: ! Private functions 
 38: PRIVATE :: PERIODIC_RESET, OPP_PER, OPP_PER_TRI, CALC_ENERGY 
 39: PRIVATE :: CALC_GRAD, CALC_FCTS, CONSTRAIN_PARAMETERS 
 40: PRIVATE :: BOX_DERIV, CALC_BOX_DERIV, CHECK_ANGLES, REJECT, CONSTRAIN_VOLUME 
 41:  
 42: ! K: frequency parameter, controls the number of wells 
 43: ! PHI: phase parameter, controls the position of the first minimum 
 44: ! RCUT: cutoff distance at which the energy, gradients and Hessian of 
 45: !       the potential go smoothly to zero. 
 46: DOUBLE PRECISION :: OPP_K, OPP_PHI, OPP_RCUT 
 47:  
 48: ! OPP_MODE: 0 standard minimisation 
 49: !           1 the final triplet of coordinates are the unit cell length 
 50: !             parameters and the second last triplet are the unit cell 
 51: !             angles, giving a triclinic unit cell, which will be optimised 
 52: !           2 as for 1, but the third last triplet are the potential 
 53: !             parameters OPP_K, OPP_PHI and 0 and 
 54: !             the potential parameters will be optimised, depending on the 
 55: !             value of OPP_PARAMS 
 56: ! OPP_PARAMS: Specifies which of the three potential parameters will be 
 57: !             optimised. An integer between 1 and 3. 
 58: !           1 OPP_K will be optimised 
 59: !           2 OPP_PHI will be optimised 
 60: !           3 both parameters will be optimised 
 61: INTEGER :: OPP_MODE, OPP_PARAMS 
 62:  
 63: ! V_RCUT: the unmodified potential evaluated at the cutoff 
 64: ! FD_RCUT: the first derivative evaluated at the cutoff 
 65: ! SD_RCUT: the second derivative evaluated at the cutoff 
 66: DOUBLE PRECISION :: V_RCUT, FD_RCUT!, SD_RCUT 
 67:  
 68: ! PI: the usual meaning 
 69: DOUBLE PRECISION, PARAMETER :: PI = 4*ATAN(1.D0) 
 70:  
 71: ! IMAGE_CUTOFF: if more periodic images than this in any direction fall 
 72: !               within the cutoff, the evaluation will not be carried 
 73: !               out and the step will be rejected. 
 74: DOUBLE PRECISION, PARAMETER :: IMAGE_CUTOFF = 30 
 75:  
 76: ! Parameters specifying the largest steps in OPP_TAKESTEP 
 77: ! MAX_LENGTH_STEP: largest allowed step in the length of a lattice 
 78: !                  vector 
 79: ! MAX_ANGLE_STEP: largest allowed step in the angle of a lattice vector 
 80: ! MAX_K_STEP: largest allowed step in OPP_K 
 81: ! MAX_PHI_STEP: largest allowed step in OPP_PHI 
 82: DOUBLE PRECISION, PARAMETER :: MAX_LENGTH_STEP = 0.3D0 
 83: DOUBLE PRECISION, PARAMETER :: MAX_ANGLE_STEP = 0.1D0 
 84: DOUBLE PRECISION, PARAMETER :: MAX_K_STEP = 1.D0 
 85: DOUBLE PRECISION, PARAMETER :: MAX_PHI_STEP = 0.5D0 
 86:  
 87: ! Parameters used for constraining the unit cell volume and the potential 
 88: ! parameters 
 89: ! V_EPS: scaling factor for the unit cell volume contraint energy 
 90: ! V_SIGMA: WCA cutoff distance 
 91: ! K_LOWER, K_UPPER: bounds on OPP_K 
 92: ! K_REPULSION: harmonic force constant for OPP_K 
 93: ! OPP_PHI is periodic over the range 0 to 2*PI, so doesn't need constraints 
 94: DOUBLE PRECISION, PARAMETER :: V_EPS = 1.D-3 
 95: DOUBLE PRECISION, PARAMETER :: V_SIGMA = 3.D-1 
 96: DOUBLE PRECISION, PARAMETER :: K_LOWER = 5.D0, K_UPPER = 20.D0 
 97: DOUBLE PRECISION, PARAMETER :: K_REPULSION = 1.D4 
 98:  
 99: ! BOX_DERIV: stores information for calculating lattice derivatives 
100: TYPE :: BOX_DERIV 
101:     ! ORTHOG: the orthogonal tansformation matrix, for changing from 
102:     !         crystallographic space to orthogonal space first index is row, 
103:     !         second is column 
104:     DOUBLE PRECISION :: ORTHOG(3,3) 
105:  
106:     ! DERIV: derivatives of the orthogonal matrix wrt lattice parameters, 
107:     !        first index indicates the row, second index the column and 
108:     !        third index the parameter, with 1-3 being the angles and 4-6 being 
109:     !        the lengths 
110:     DOUBLE PRECISION :: DERIV(3, 3, 6) 
111:  
112:     ! RECIP: lengths of the reciprocal lattice vectors 
113:     DOUBLE PRECISION :: RECIP(3) 
114:  
115:     ! CELL_RANGE: number of cells in each direction that must be evaluated 
116:     INTEGER :: CELL_RANGE(3) 
117:  
118:     ! V: dimensionless volume of the unit cell, ie volume if all lengths were 1 
119:     DOUBLE PRECISION :: V 
120:  
121:     ! V_DERIV: derivatievs of V wrt the lattice angles 
122:     DOUBLE PRECISION :: V_DERIV(3) 
123: END TYPE BOX_DERIV 
124:  
125: CONTAINS 
126:  
127: !------------------------------------------------------------------------------- 
128: ! 
129: ! Main routine for calculating the energy and gradients. 
130: ! X: coordinates array 
131: ! GRAD: gradients array 
132: ! ENERGY: calculated energy 
133: ! GTEST: whether gradients should be calculated 
134: ! 
135: !------------------------------------------------------------------------------- 
136:     SUBROUTINE OPP(X, GRAD, ENERGY, GTEST) 
137:  
138:         ! MYUNIT: file unit for main output file 
139:         ! NATOMS: number of particles 
140:         ! PERIODIC: whether to use periodic boundary conditions 
141:         USE COMMONS, ONLY: MYUNIT, NATOMS, PERIODIC 
142:  
143:         IMPLICIT NONE 
144:  
145:         DOUBLE PRECISION, INTENT(INOUT) :: X(3*NATOMS) 
146:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOMS), ENERGY 
147:         LOGICAL, INTENT(IN) :: GTEST 
148:  
149:         ! I, J: indices for looping over particles 
150:         ! NMOL: actual number of particles 
151:         INTEGER :: I, J, NMOL 
152:  
153:         ! RIJ: Interparticle vector 
154:         ! MODRIJ: distance between a pair of particles 
155:         ! DVDR: derivative of pair potential wrt MODRIJ 
156:         DOUBLE PRECISION :: RIJ(3), MODRIJ, DVDR 
157:  
158:         ! Factors for triclinic lattice 
159:         TYPE(BOX_DERIV) :: BD 
160:  
161:         ! Initialise output variables 
162:         ENERGY = 0.D0 
163:         GRAD(:) = 0.D0 
164:  
165:         ! Calculate factors that do depend on parameters, but not on RIJ 
166:         IF (OPP_MODE .EQ. 2) THEN 
167:             CALL CALC_FCTS(X(NATOMS*3-8:NATOMS*3-6)) 
168:         ELSE 
169:             CALL CALC_FCTS() 
170:         END IF 
171:  
172:         IF (.NOT. PERIODIC) THEN 
173:             IF (OPP_MODE .EQ. 0) THEN 
174:                 ! Normal cluster 
175:                 NMOL = NATOMS 
176:                 DO I = 1, NMOL-1 ! Outer loop over particles 
177:                     DO J = I+1, NMOL ! Inner loop over particles 
178:                         ! Get the particle-particle distance 
179:                         RIJ(:) = X(3*I-2:3*I) - X(3*J-2:3*J) 
180:                         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:))) 
181:  
182:                         ! Add energy and gradients 
183:                         IF (MODRIJ .LT. OPP_RCUT) THEN 
184:                             ENERGY = ENERGY + CALC_ENERGY(MODRIJ) 
185:  
186:                             IF (GTEST) THEN 
187:                                 DVDR = CALC_GRAD(MODRIJ) 
188:                                 GRAD(3*I-2:3*I) = GRAD(3*I-2:3*I) + & 
189:                                                   DVDR*RIJ(:)/MODRIJ 
190:                                 GRAD(3*J-2:3*J) = GRAD(3*J-2:3*J) - & 
191:                                                   DVDR*RIJ(:)/MODRIJ 
192:                             END IF ! GTEST 
193:                         END IF ! If less than cutoff distance 
194:                     END DO ! Inner loop over particles 
195:                 END DO ! Outer loop over particles 
196:             ELSE ! mode not equal to zero 
197:                 WRITE (MYUNIT, *) 'OPP> ERROR, PERIODIC must be ', & 
198:                                   'specified if mode is not 0' 
199:                 STOP 
200:             END IF 
201:         ELSE ! PERIODIC 
202:             SELECT CASE (OPP_MODE) 
203:             CASE(0) 
204:                 ! Periodic, fixed box 
205:                 ! Reset all particles to within the box 
206:                 CALL PERIODIC_RESET(X(:)) 
207:  
208:                 NMOL = NATOMS 
209:                 ! We need to include self-interactions of particles in different 
210:                 ! unit cells 
211:                 DO I = 1, NMOL ! Outer loop over particles 
212:                     DO J = I, NMOL ! Inner loop over particles 
213:                         ! Add energy and gradients 
214:                         CALL OPP_PER(I, J, X(3*I-2:3*I), X(3*J-2:3*J), & 
215:                                      GRAD(3*I-2:3*I), GRAD(3*J-2:3*J), & 
216:                                      ENERGY, GTEST) 
217:                     END DO ! Inner loop over particles 
218:                 END DO ! Outer loop over particles 
219:  
220:             CASE(1) 
221:                 ! Triclinic varying lattice 
222:                 ! Reset all particles to within the box 
223:                 CALL PERIODIC_RESET(X(:)) 
224:  
225: !                DO I = 1, NATOMS 
226: !                    WRITE (MYUNIT, *) X(3*I-2:3*I) 
227: !                END DO 
228:  
229:                 ! Calculate box derivative factors 
230:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), & 
231:                                     X(3*NATOMS-2:3*NATOMS  )) 
232:  
233:                 ! Check whether the unit cell angles are physically possible. 
234:                 ! Reject if not. 
235:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN 
236:                     CALL REJECT(ENERGY, GRAD) 
237:                     RETURN 
238:                 END IF 
239:  
240:                 ! Reject steps where the cell range has got bigger than the  
241:                 ! cutoff. Such unit cells are highly distorted and there are  
242:                 ! probably equivalent ones of a more usual shape. Allowing them 
243:                 ! leads to very slow run times. 
244:                 IF (.NOT. ALL(BD%CELL_RANGE .LE. IMAGE_CUTOFF)) THEN 
245:                     CALL REJECT(ENERGY, GRAD) 
246:                     RETURN 
247:                 END IF 
248:  
249:                 NMOL = NATOMS-2 
250:                 ! We need to include self-interactions of particles in different 
251:                 ! unit cells 
252:                 DO I = 1, NMOL! Outer loop over particles 
253:                     DO J = I, NMOL ! Inner loop over particles 
254:                         ! Add energy and gradients 
255:                         CALL OPP_PER_TRI(I, J, X(:), GRAD(3*I-2:3*I), & 
256:                             GRAD(3*J-2:3*J), ENERGY, GTEST, & 
257:                             GRAD(3*NATOMS-5:3*NATOMS), BD) 
258:                     END DO ! Inner loop over particles 
259:                 END DO ! Outer loop over particles 
260:  
261:                 ! Constrain the box volume to be greater than 0 with a WCA-style 
262:                 ! repulsion 
263:                 CALL CONSTRAIN_VOLUME(ENERGY, GRAD(3*NATOMS-5:3*NATOMS-3), BD, & 
264:                                       GTEST) 
265:  
266:             CASE(2) 
267:                 ! Triclinic varying lattice, with varying parameters 
268:                 ! Reset all particles to within the box 
269:                 CALL PERIODIC_RESET(X(:)) 
270:  
271:                 ! Calculate box derivative factors 
272:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), & 
273:                                     X(3*NATOMS-2:3*NATOMS  )) 
274:  
275:                 ! Check whether the unit cell angles are physically possible. 
276:                 ! If we have a disallowed unit cell, we set the energy to very 
277:                 ! high and the gradients to very small, so the step is 
278:                 ! 'converged' but gets rejected. 
279:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN 
280:                     CALL REJECT(ENERGY, GRAD) 
281:                 END IF 
282:  
283:                 ! Reject steps where the cell range has got bigger than the  
284:                 ! cutoff. Such unit cells are highly distorted and there are  
285:                 ! probably equivalent ones of a more usual shape. Allowing them 
286:                 ! leads to very slow run times. 
287:                 IF (.NOT. ALL(BD%CELL_RANGE .LE. IMAGE_CUTOFF)) THEN 
288:                     CALL REJECT(ENERGY, GRAD) 
289:                     RETURN 
290:                 END IF 
291:  
292:                 NMOL = NATOMS-3 
293:                 ! We need to include self-interactions of particles in different 
294:                 ! unit cells 
295:                 DO I = 1, NMOL! Outer loop over particles 
296:                     DO J = I, NMOL ! Inner loop over particles 
297:                         ! Add energy and gradients 
298:                         CALL OPP_PER_TRI(I, J, X(:), GRAD(3*I-2:3*I), & 
299:                             GRAD(3*J-2:3*J), ENERGY, GTEST, & 
300:                             GRAD(3*NATOMS-5:3*NATOMS), BD, & 
301:                             GRAD(3*NATOMS-8:3*NATOMS-6)) 
302:                     END DO ! Inner loop over particles 
303:                 END DO ! Outer loop over particles 
304:  
305:                 ! Constrain the box volume to be greater than 0 with a WCA-style 
306:                 ! repulsion 
307:                 CALL CONSTRAIN_VOLUME(ENERGY, GRAD(3*NATOMS-5:3*NATOMS-3), BD, & 
308:                                       GTEST) 
309:  
310:                 ! Constrain the potential parameters with a harmonic repulsion 
311:                 CALL CONSTRAIN_PARAMETERS(ENERGY, GRAD(3*NATOMS-8:3*NATOMS-6), & 
312:                                           GTEST) 
313:             END SELECT ! Mode selections 
314:         END IF ! Periodic 
315:  
316:         FLUSH(MYUNIT) 
317:  
318:     END SUBROUTINE OPP 
319:  
320: !------------------------------------------------------------------------------- 
321: ! 
322: ! Rejects a step by setting the energy very high and the gradients very low. 
323: ! This tricks L-BFGS into thinking it's converged, but the high energy of the 
324: ! quench will ensure it is rejected. 
325: ! ENERGY: system energy 
326: ! GRAD: system gradients 
327: ! 
328: !------------------------------------------------------------------------------- 
329:  
330:     SUBROUTINE REJECT(ENERGY, GRAD) 
331:  
332:         ! NATOMS: number of atoms 
333:         USE COMMONS, ONLY: NATOMS 
334:  
335:         IMPLICIT NONE 
336:  
337:         DOUBLE PRECISION, INTENT(OUT) :: ENERGY, GRAD(3*NATOMS) 
338:  
339:         ENERGY = 1.D20 
340:         GRAD(:) = 1.D-20 
341:  
342:     END SUBROUTINE REJECT 
343:  
344: !------------------------------------------------------------------------------- 
345: ! 
346: ! Calculate factors that do depend on parameters, but not on RIJ 
347: ! Some book keeping of the potential parameters to deal with different modes 
348: ! PARAMS: the potential parameters when those are to be adjusted, otherwise 
349: !         ignored 
350: ! 
351: !------------------------------------------------------------------------------- 
352:  
353:     SUBROUTINE CALC_FCTS(PARAMS) 
354:  
355:         IMPLICIT NONE 
356:  
357:         DOUBLE PRECISION, INTENT(INOUT), OPTIONAL :: PARAMS(3) 
358:  
359:         ! S, C: factors which occur repeatedly 
360:         DOUBLE PRECISION :: S, C 
361:  
362:         ! Adjust the potential parameters, if the mode requires it 
363:         IF (OPP_MODE .EQ. 2 .AND. PRESENT(PARAMS)) THEN 
364:             ! For each parameter, if it is varying, set it from the coordinate. 
365:             ! If it's not varying, reset the coordinate, which TAKESTEP will 
366:             ! have changed. 
367:             IF (BTEST(OPP_PARAMS, 0)) THEN 
368:                 OPP_K = PARAMS(1) 
369:             ELSE 
370:                 PARAMS(1) = OPP_K 
371:             END IF 
372:             IF (BTEST(OPP_PARAMS, 1)) THEN 
373:                 OPP_PHI = PARAMS(2) 
374:             ELSE 
375:                 PARAMS(2) = OPP_PHI 
376:             END IF 
377:             PARAMS(3) = 0.D0 
378:         END IF 
379:  
380:         ! Calculate factors 
381:         S = SIN(OPP_K * (OPP_RCUT - 1) - OPP_PHI) 
382:         C = COS(OPP_K * (OPP_RCUT - 1) - OPP_PHI) 
383:  
384:         ! Unmodifed potential evaluated at the cutoff 
385:         V_RCUT = OPP_RCUT**(-15) + OPP_RCUT**(-3) * C 
386:         ! First derivative evaluated at the cutoff 
387:         FD_RCUT = - 15.D0 * OPP_RCUT**(-16) - OPP_K * OPP_RCUT**(-3) * S - & 
388:                     3.D0 * OPP_RCUT**(-4) * C 
389:         ! Second derivative evaluated at the cutoff 
390: !        SD_RCUT = 240.D0 * OPP_RCUT**(-17) - OPP_K**2 * OPP_RCUT**(-3) * C + & 
391: !                  6.D0 * OPP_K * OPP_RCUT**(-4) * S + 12.D0 * OPP_RCUT**(-5) * C 
392:  
393:     END SUBROUTINE CALC_FCTS 
394:  
395: !------------------------------------------------------------------------------- 
396: ! 
397: ! Calculates the energy contribution for a pair of particles. 
398: ! MODRIJ: distance beteen the particle pair 
399: ! 
400: !------------------------------------------------------------------------------- 
401:     PURE DOUBLE PRECISION FUNCTION CALC_ENERGY(MODRIJ) 
402:  
403:         IMPLICIT NONE 
404:  
405:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
406:  
407:         ! ENERGY: working variable for energy 
408:         DOUBLE PRECISION :: ENERGY 
409:  
410:         ! Unmodified potential 
411:         ENERGY = MODRIJ**(-15) + MODRIJ**(-3) * & 
412:                  COS(OPP_K * (MODRIJ - 1) - OPP_PHI) 
413:  
414:         ! Term to make the potential go to zero at OPP_RCUT 
415:         ENERGY = ENERGY - V_RCUT 
416:  
417:         ! Term to make the first derivative go to zero at OPP_RCUT 
418:         ENERGY = ENERGY - (MODRIJ - OPP_RCUT) * FD_RCUT 
419:  
420:         ! Term to make the second derivative go to zero at OPP_RCUT 
421: !        ENERGY = ENERGY - 0.5D0 * (MODRIJ - OPP_RCUT)**2 * SD_RCUT 
422:  
423:         ! Set the output 
424:         CALC_ENERGY = ENERGY 
425:  
426:     END FUNCTION CALC_ENERGY 
427:  
428: !------------------------------------------------------------------------------- 
429: ! 
430: ! Calculates the pair potential gradient for a pair of particles. 
431: ! MODRIJ: distance beteen the particle pair 
432: ! 
433: !------------------------------------------------------------------------------- 
434:     PURE DOUBLE PRECISION FUNCTION CALC_GRAD (MODRIJ) 
435:  
436:         IMPLICIT NONE 
437:  
438:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
439:  
440:         ! DVDR: temporary store for the result 
441:         DOUBLE PRECISION :: DVDR 
442:  
443:         ! Unmodified gradient 
444:         DVDR = - 15.D0 * MODRIJ**(-16) - & 
445:                OPP_K * MODRIJ**(-3) * SIN(OPP_K * (MODRIJ - 1) - OPP_PHI) - &  
446:                3.D0 * MODRIJ**(-4) * COS(OPP_K * (MODRIJ - 1) - OPP_PHI) 
447:  
448:         ! Term to make the first derivative go to zero at OPP_RCUT 
449:         DVDR = DVDR - FD_RCUT 
450:  
451:         ! Term to make the second derivative go to zero at OPP_RCUT 
452:         DVDR = DVDR - (MODRIJ - OPP_RCUT) * SD_RCUT 
453:  
454:         CALC_GRAD = DVDR 
455:  
456:     END FUNCTION CALC_GRAD 
457:  
458: !------------------------------------------------------------------------------- 
459: ! 
460: ! Initialisation. Check values for OPP_MODE and OPP_PARAMS. 
461: ! 
462: !------------------------------------------------------------------------------- 
463:     SUBROUTINE OPP_INITIALISE() 
464:  
465:         ! MYUNIT: file unit for main GMIN output 
466:         USE COMMONS, ONLY: MYUNIT 
467:  
468:         IMPLICIT NONE 
469:  
470:         ! Check we have sensible value for the mode 
471:         IF (OPP_MODE .NE. 0 .AND. OPP_MODE .NE. 1 .AND. & 
472:             OPP_MODE .NE. 2 ) THEN 
473:             WRITE (MYUNIT, *) 'OPP> ERROR: mode must be 0, 1, 2' 
474:             STOP 
475:         END IF 
476:  
477:         ! Check we have sensible values for the params 
478:         IF (OPP_MODE .EQ. 2) THEN 
479:             IF (OPP_PARAMS .LT. 1 .OR. OPP_PARAMS .GT. 3) THEN 
480:                 WRITE (MYUNIT, *) 'OPP> ERROR: params must be between 1 and 3' 
481:                 STOP 
482:             END IF 
483:         END IF           
484:  
485:     END SUBROUTINE OPP_INITIALISE 
486:  
487: !------------------------------------------------------------------------------- 
488: ! 
489: ! Takes a step, perturbing the Cartesian coordinates, the unit cell parameters 
490: ! and the potential parameters, according to the specified mode. 
491: ! NP: the index in the main COORDS array to take the step from 
492: ! 
493: !------------------------------------------------------------------------------- 
494:  
495:     SUBROUTINE OPP_TAKESTEP(NP) 
496:  
497:         ! COORDS: full set of Markov chain coordinates 
498:         ! MYUNIT: file unit of the main GMIN output file 
499:         ! NATOMS: number of coordinates 
500:         ! PERIODIC: whetehr periodic boundary conditions are used 
501:         ! PERCOLATET: whether percolation is used to constrain particles 
502:         ! RADIUS: container radius 
503:         ! STEP: maximum step size for this move 
504:         ! TMOVE: specifies which step to do a translational move on 
505:         ! TWOD: whetehr the system is two dimensional 
506:         USE COMMONS, ONLY: COORDS, MYUNIT, NATOMS, PERIODIC, PERCOLATET, &  
507:                            RADIUS, STEP, TMOVE, TWOD 
508:          
509:         ! VEC_LEN: function, returns length of a 3D vector 
510:         ! VEC_RANDOM: function, returns a randomly orientated 3D unit vector 
511:         USE VEC3, ONLY: VEC_LEN, VEC_RANDOM 
512:  
513:         IMPLICIT NONE 
514:  
515:         INTEGER, INTENT(IN) :: NP 
516:  
517:         ! X: local copy of the coordinates 
518:         ! RAND_STEP: unit vector in a random step direction 
519:         ! STEP_SIZE: random fraction of the max step size to move 
520:         ! NE_ANGLES: new set of unit cell angles, to be checked 
521:         DOUBLE PRECISION :: X(3*NATOMS), RAND_STEP(3), STEP_SIZE, NEW_ANGLES(3) 
522:  
523:         ! I: iteration index 
524:         ! NMOL: actual number of atoms 
525:         INTEGER :: I, NMOL 
526:  
527:         ! Set the local copy of the coordinates 
528:         X(:) = COORDS(:, NP) 
529:  
530:         ! Set NMOL 
531:         NMOL = NATOMS - OPP_MODE - 1 
532:  
533:         ! Random translational steps for each of the particles 
534:         IF (TMOVE(NP)) THEN 
535:             DO I = 1, NMOL 
536:                 ! Generate either a 2D or 3D random vector 
537:                 ! The magnitude will be uniformly distributed in a circle or  
538:                 ! sphere with radius given by the value of STEP for this step 
539:                 IF (TWOD) THEN 
540:                     DO 
541:                         CALL RANDOM_NUMBER(RAND_STEP(1)) 
542:                         CALL RANDOM_NUMBER(RAND_STEP(2)) 
543:                         RAND_STEP(3) = 0.D0 
544:                         IF(VEC_LEN(RAND_STEP) .LE. 1.D0) EXIT 
545:                     END DO 
546:                     RAND_STEP(:) = RAND_STEP(:) * STEP(NP) 
547:                 ELSE ! 3D 
548:                     RAND_STEP(:) = VEC_RANDOM() 
549:                     CALL RANDOM_NUMBER(STEP_SIZE) 
550:                     ! Sqrt for random volume distribution 
551:                     RAND_STEP(:) = RAND_STEP(:) * SQRT(STEP_SIZE * STEP(NP)) 
552:                 END IF 
553:                 ! Add the geenrated ste onto the coordinates 
554:                 X(I*3-2:I*3) = X(I*3-2:I*3) + RAND_STEP(:) 
555:             END DO 
556:         END IF 
557:  
558:         ! If not periodic and we're using a radius container, bring particles in 
559:         IF (.NOT. PERIODIC .AND. .NOT. PERCOLATET) THEN 
560:             DO I = 1, NMOL 
561:                 IF (VEC_LEN(X(I*3-2:I*3)) .GT. RADIUS) THEN 
562:                     WRITE(MYUNIT, *) 'OPP_TAKESTEP> ', & 
563:                         'coord outside container, bringing in' 
564:                     X(I*3-2:I*3) = X(I*3-2:I*3) - & 
565:                         SQRT(RADIUS) * NINT(X(I*3-2:I*3) / SQRT(RADIUS)) 
566:                 END IF 
567:             END DO 
568:         END IF 
569:  
570:         ! Box length steps 
571:         IF (OPP_MODE .GE. 1) THEN 
572:             X(3*NATOMS-2:3*NATOMS) = X(3*NATOMS-2:3*NATOMS) + & 
573:                                      VEC_RANDOM() * MAX_LENGTH_STEP 
574:         END IF 
575:  
576:         ! Box angle steps 
577:         IF (OPP_MODE .GE. 1) THEN 
578:             DO 
579:                 NEW_ANGLES(:) = X(3*NATOMS-5:3*NATOMS-3) + & 
580:                                 VEC_RANDOM() * MAX_ANGLE_STEP 
581:  
582:                 ! See whether the unit cell we've generated is valid 
583:                 IF (CHECK_ANGLES(NEW_ANGLES(:))) THEN 
584:                     X(3*NATOMS-5:3*NATOMS-3) =  NEW_ANGLES(:) 
585:                     EXIT 
586:                 END IF 
587:             END DO 
588:         END IF 
589:  
590:         ! Parameter steps 
591:         IF (OPP_MODE .EQ. 2) THEN 
592:             IF (BTEST(OPP_PARAMS, 0)) THEN 
593:                 CALL RANDOM_NUMBER(STEP_SIZE) 
594:                 X(3*NATOMS-8) = X(3*NATOMS-8) + & 
595:                                 (STEP_SIZE - 0.5D0) * 2.D0 * MAX_K_STEP 
596:             END IF 
597:             IF (BTEST(OPP_PARAMS, 1)) THEN 
598:                 CALL RANDOM_NUMBER(STEP_SIZE) 
599:                 X(3*NATOMS-7) = X(3*NATOMS-7) + & 
600:                                 (STEP_SIZE - 0.5D0) * 2.D0 * MAX_PHI_STEP 
601:             END IF             
602:         END IF 
603:  
604:         COORDS(:, NP) = X(:) 
605:  
606:     END SUBROUTINE OPP_TAKESTEP 
607:  
608: !------------------------------------------------------------------------------- 
609: ! 
610: ! Resets all coordinates to within the unit cell 
611: ! The unit cell runs from 0 to 1 in all directions. Coordinates are later scaled 
612: ! by the box lengths. 
613: ! COORDS: coordinates of the particles 
614: ! 
615: !------------------------------------------------------------------------------- 
616:     SUBROUTINE PERIODIC_RESET(COORDS) 
617:  
618:         ! BOXLX, BOXLY, BOXLZ: dimensions of box 
619:         ! NATOMS: number of particles 
620:         USE COMMONS, ONLY: BOXLX, BOXLY, BOXLZ, NATOMS 
621:  
622:         IMPLICIT NONE 
623:  
624:         DOUBLE PRECISION, INTENT(INOUT) :: COORDS(3*NATOMS) 
625:  
626:         ! Sort out the unit cell sizes, if they are varying 
627:         IF (OPP_MODE .GE. 1) THEN 
628:             BOXLX = COORDS(3*NATOMS - 2) 
629:             BOXLY = COORDS(3*NATOMS - 1) 
630:             BOXLZ = COORDS(3*NATOMS) 
631:         END IF 
632:  
633:         ! Reset coordinates to within the box 
634:         SELECT CASE (OPP_MODE) 
635:         CASE(0) 
636:             COORDS(:) = COORDS(:) - FLOOR(COORDS(:)) 
637:         CASE(1) 
638:             COORDS(1:3*NATOMS-6) = COORDS(1:3*NATOMS-6) - & 
639:                                    FLOOR(COORDS(1:3*NATOMS-6)) 
640:             ! Make the lattice angles modulo 2*pi 
641:             ! They will be checked later 
642:             COORDS(3*NATOMS-5:3*NATOMS-3) = COORDS(3*NATOMS-5:3*NATOMS-3) - & 
643:                 2.D0*PI*FLOOR(COORDS(3*NATOMS-5:3*NATOMS-3)/(2.D0*PI)) 
644:         CASE(2) 
645:            COORDS(1:3*NATOMS-9) = COORDS(1:3*NATOMS-9) - & 
646:                                    FLOOR(COORDS(1:3*NATOMS-9)) 
647:             ! Make the lattice angles modulo 2*pi 
648:             ! They will be checked later 
649:             COORDS(3*NATOMS-5:3*NATOMS-3) = COORDS(3*NATOMS-5:3*NATOMS-3) - & 
650:                 2.D0*PI*FLOOR(COORDS(3*NATOMS-5:3*NATOMS-3)/(2.D0*PI)) 
651:             ! Make OPP_PHI modulo 2*pi 
652:             COORDS(3*NATOMS-7) = COORDS(3*NATOMS-7) - & 
653:                 2.D0*PI*FLOOR(COORDS(3*NATOMS-7)/(2.D0*PI)) 
654:         END SELECT 
655:  
656:     END SUBROUTINE PERIODIC_RESET 
657:  
658: !------------------------------------------------------------------------------- 
659: ! 
660: ! Finds the energy and gradients for the othogonal fixed periodic system 
661: ! IDI, IDJ: id numbers of the particles 
662: ! COORDS_I: coordinates of particle I 
663: ! COORDS_J: coordinates of particle J 
664: ! GRAD_I: gradients for particle I 
665: ! GRAD_J: gradients for particle J 
666: ! ENERGY: energy contribution for this pair 
667: ! GTEST: whether to calculate gradients 
668: ! 
669: !------------------------------------------------------------------------------- 
670:  
671:     SUBROUTINE OPP_PER(IDI, IDJ, COORDS_I,COORDS_J,GRAD_I,GRAD_J,ENERGY, GTEST) 
672:  
673:         ! BOXLX, BOXLY, BOXLZ: dimensions of the box 
674:         USE COMMONS, ONLY: BOXLX, BOXLY, BOXLZ 
675:  
676:         IMPLICIT NONE 
677:  
678:         INTEGER, INTENT(IN) :: IDI, IDJ 
679:         DOUBLE PRECISION, INTENT(IN) :: COORDS_I(3), COORDS_J(3) 
680:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_I(3), GRAD_J(3), ENERGY 
681:         LOGICAL, INTENT(IN) :: GTEST 
682:  
683:         ! XJ: copy of j coordinates 
684:         ! RIJ: particle-particle vector 
685:         ! MODRIJ: particle-particle distance 
686:         ! DVDR: derivative of pair potential wrt MODRIJ 
687:         ! BOX_SIZE: convenience copy of the box size 
688:         DOUBLE PRECISION :: XJ(3), RIJ(3), MODRIJ, DVDR, BOX_SIZE(3) 
689:  
690:         ! PER_k: integers for looping over cell possibilities 
691:         ! CELL_RANGE: number of cells in each direction required 
692:         INTEGER :: PER_X, PER_Y, PER_Z, CELL_RANGE(3) 
693:  
694: !        DOUBLE PRECISION :: START_ENERGY 
695: !        START_ENERGY = ENERGY 
696:  
697:         BOX_SIZE(1) = BOXLX 
698:         BOX_SIZE(2) = BOXLY 
699:         BOX_SIZE(3) = BOXLZ 
700:  
701:         ! Set the number of cells we need to check, based on the cell size and 
702:         ! the potential cutoff. 
703:         CELL_RANGE(:) = CEILING(OPP_RCUT/BOX_SIZE(:)) 
704:  
705:         ! Loop over different coordinate possibilities. 
706:         DO PER_X = - CELL_RANGE(1), CELL_RANGE(1) 
707:             DO PER_Y = - CELL_RANGE(2), CELL_RANGE(2) 
708:                 DO PER_Z = - CELL_RANGE(3), CELL_RANGE(3) 
709:  
710:                     ! Skip the interaction of a particle with itself in the same 
711:                     ! cell and don't double count interactions within the unit 
712:                     ! cell. 
713:                     IF (PER_X .EQ. 0 .AND. PER_Y .EQ. 0 .AND. PER_Z .EQ. 0 & 
714:                         .AND. IDI .EQ. IDJ) CYCLE 
715:  
716:                     ! Adjust the j coordinates for the periodic image 
717:                     XJ(1) = COORDS_J(1) + PER_X 
718:                     XJ(2) = COORDS_J(2) + PER_Y 
719:                     XJ(3) = COORDS_J(3) + PER_Z 
720:  
721:                     ! Find the distance and scale by the box size 
722:                     RIJ(:) = (COORDS_I(:) - XJ(:))*BOX_SIZE(:) 
723:                     MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:))) 
724:  
725:                     IF (MODRIJ .LT. OPP_RCUT) THEN 
726:                         ! Add energy and gradients 
727:                         ! Divide by 2 for images of the same atom, to avoid 
728:                         ! double counting 
729:                         ENERGY = ENERGY + MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * & 
730:                                           CALC_ENERGY(MODRIJ) 
731:  
732:                         IF (GTEST) THEN 
733:                             ! Divide gradient by 2 for images of the same atom, 
734:                             ! to avoid double counting 
735:                             DVDR = MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * & 
736:                                    CALC_GRAD(MODRIJ) 
737:                             ! We must undo the box size coordinate scaling 
738:                             GRAD_I(:) = GRAD_I(:) + & 
739:                                         DVDR*RIJ(:)*BOX_SIZE(:)/MODRIJ 
740:                             GRAD_J(:) = GRAD_J(:) - & 
741:                                         DVDR*RIJ(:)*BOX_SIZE(:)/MODRIJ 
742:                         END IF ! GTEST 
743:                     END IF ! End if less than cutoff 
744:                 END DO ! z loop 
745:             END DO ! y loop 
746:         END DO ! x loop 
747:  
748:     END SUBROUTINE OPP_PER 
749:  
750: !------------------------------------------------------------------------------- 
751: ! 
752: ! Calculates factors for the lattice derivatives that only depend on the cell 
753: ! parameters. 
754: ! ANGLES: alpha, beta and gamma lattice parameters 
755: ! SIZES: a, b and c lattice parameters 
756: ! 
757: !------------------------------------------------------------------------------- 
758:  
759:     PURE TYPE(BOX_DERIV) FUNCTION CALC_BOX_DERIV(ANGLES, SIZES) 
760:  
761:         IMPLICIT NONE 
762:  
763:         DOUBLE PRECISION, INTENT(IN) :: ANGLES(3), SIZES(3) 
764:  
765:         ! V: factor, related to (but not quite equal to) the unit cell volume 
766:         ! C, S: cos and sin of the lattice angles 
767:         ! ORTHOG, DERIV: parts of the output 
768:         DOUBLE PRECISION :: V, C(3), S(3) 
769:         DOUBLE PRECISION :: ORTHOG(3, 3), DERIV(3, 3, 6) 
770:  
771:         ! Initialise output variables 
772:         ORTHOG(:,:) = 0.D0 
773:         DERIV(:,:,:) = 0.D0 
774:  
775:         ! Calculate factors 
776:         C(:) = COS(ANGLES(:)) 
777:         S(:) = SIN(ANGLES(:)) 
778:         V = SQRT(1 - C(1)**2 - C(2)**2 - C(3)**2 + 2.D0 * C(1) * C(2) * C(3)) 
779:  
780:         ! Calculate the orthogonal transformation matrix 
781:         ORTHOG(1, 1) = SIZES(1) 
782:         ORTHOG(1, 2) = SIZES(2) * C(3) 
783:         ORTHOG(1, 3) = SIZES(3) * C(2) 
784:         ORTHOG(2, 2) = SIZES(2) * S(3) 
785:         ORTHOG(2, 3) = SIZES(3) * (C(1) - C(2) * C(3)) / S(3) 
786:         ORTHOG(3, 3) = SIZES(3) * V / S(3) 
787:  
788:         ! Calculate the derivatives of the othogonal matrix 
789:         ! Derivatives of the top row wrt angles 
790:         DERIV(1, 3, 2) = - SIZES(3) * S(2) 
791:         DERIV(1, 2, 3) = - SIZES(2) * S(3) 
792:         ! Derivatives of the top row wrt to lengths 
793:         DERIV(1, 1, 4) = 1.D0 
794:         DERIV(1, 2, 5) = C(3) 
795:         DERIV(1, 3, 6) = C(2) 
796:         ! Derivatives of the middle row wrt angles 
797:         DERIV(2, 3, 1) = - SIZES(3) * S(1) / S(3) 
798:         DERIV(2, 3, 2) = SIZES(3) * S(2) * C(3) / S(3) 
799:         DERIV(2, 2, 3) = SIZES(2) * C(3) 
800:         DERIV(2, 3, 3) = SIZES(3) * (C(2) - C(1) * C(3)) / S(3)**2 
801:         ! Derivatives of the middle row wrt to lengths 
802:         DERIV(2, 2, 5) = S(3) 
803:         DERIV(2, 3, 6) = (C(1) - C(2) * C(3)) / S(3) 
804:         ! Derivatives of the bottom row wrt angles 
805:         DERIV(3, 3, 1) = SIZES(3) * S(1) * (C(1) - C(2) * C(3)) / (S(3) * V) 
806:         DERIV(3, 3, 2) = SIZES(3) * S(2) * (C(2) - C(3) * C(1)) / (S(3) * V) 
807:         DERIV(3, 3, 3) = SIZES(3) * ((C(3) - C(1) * C(2))/V - C(3) * V/S(3)**2) 
808:         ! Derivatives of the bottom row wrt lengths 
809:         DERIV(3, 3, 6) = V / S(3) 
810:  
811:         ! Calculate the derivative of volume wrt the angles 
812:         CALC_BOX_DERIV%V_DERIV(1) = S(1) * (C(1) - C(2) * C(3)) / V 
813:         CALC_BOX_DERIV%V_DERIV(2) = S(2) * (C(2) - C(3) * C(1)) / V 
814:         CALC_BOX_DERIV%V_DERIV(3) = S(3) * (C(3) - C(1) * C(2)) / V 
815:  
816:         ! Calculate the lengths of the reciprocal lattice vectors 
817:         CALC_BOX_DERIV%RECIP(:) = S(:) / (SIZES(:) * V) 
818:  
819:         ! Calculate the number of cells to check in each direction 
820:         CALC_BOX_DERIV%CELL_RANGE(:) = CEILING(OPP_RCUT * & 
821:             CALC_BOX_DERIV%RECIP(:)) 
822:  
823:         ! Set the output 
824:         CALC_BOX_DERIV%ORTHOG(:,:) = ORTHOG(:,:) 
825:         CALC_BOX_DERIV%DERIV(:,:,:) = DERIV(:,:,:) 
826:         CALC_BOX_DERIV%V = V 
827:  
828:     END FUNCTION CALC_BOX_DERIV 
829:  
830: !------------------------------------------------------------------------------- 
831: ! 
832: ! Checks whether the unit cell angles are allowed. It is possible to generate 
833: ! sets of angles which do not correspond to a physically possible unit cell. 
834: ! This would manifest by an imaginary unit cell volume. See 
835: ! https://journals.iucr.org/a/issues/2011/01/00/au5114/au5114.pdf 
836: ! 'On the allowed values for the triclinic unit-cell angles', 
837: ! J. Foadi and G. Evans, Acta Cryst. (2011) A67, 93–95 
838: ! ANGLES: the three unit cell angles 
839: ! 
840: !------------------------------------------------------------------------------- 
841:  
842:     PURE LOGICAL FUNCTION CHECK_ANGLES(ANGLES) 
843:  
844:         IMPLICIT NONE 
845:  
846:         DOUBLE PRECISION, INTENT(IN) :: ANGLES(3) 
847:  
848:         ! SUMS: sums and differences of the angles 
849:         DOUBLE PRECISION :: SUMS(4) 
850:  
851:         ! Calculate the necessary sums 
852:         SUMS(1) =   ANGLES(1) + ANGLES(2) + ANGLES(3) 
853:         SUMS(2) = - ANGLES(1) + ANGLES(2) + ANGLES(3) 
854:         SUMS(3) =   ANGLES(1) - ANGLES(2) + ANGLES(3) 
855:         SUMS(4) =   ANGLES(1) + ANGLES(2) - ANGLES(3) 
856:  
857:         ! Check all sums 0 < SUM < 2 pi and all angles 0 < ANGLE < pi 
858:         CHECK_ANGLES = ALL(SUMS .GT. 0.D0) .AND. ALL(SUMS .LT. 2*PI) .AND. & 
859:                        ALL(ANGLES .GT. 0.D0) .AND. ALL(ANGLES .LT. PI) 
860:  
861:     END FUNCTION CHECK_ANGLES 
862:  
863: !------------------------------------------------------------------------------- 
864: ! 
865: ! Finds the energy and gradients for the periodic system with a triclinic 
866: ! varying unit cell. Optionally with varying potential parameters 
867: ! IDI, IDJ: id numbers of the particles 
868: ! COORDS: coordinates of particles (in crystallographic space) 
869: !         and lattice parameters 
870: ! GRAD_I: gradients for particle I 
871: ! GRAD_J: gradients for particle J 
872: ! ENERGY: energy contribution for this pair 
873: ! RIJ: mininum length vector between the particles 
874: ! GTEST: whether to calculate gradients 
875: ! GRAD_BOX: gradients for the box (angles, then lengths) 
876: ! BD: factors for the lattice derivatives 
877: ! GRAD_PARAMS: optional, derivatives of potential parameters 
878: ! 
879: !------------------------------------------------------------------------------- 
880:  
881:     SUBROUTINE OPP_PER_TRI(IDI, IDJ, COORDS, GRAD_I, GRAD_J, ENERGY, GTEST, & 
882:                            GRAD_BOX, BD, GRAD_PARAMS) 
883:  
884:         ! NATOMS: number of coordinates (including lattice parameters) 
885:         USE COMMONS, ONLY: NATOMS 
886:  
887:         IMPLICIT NONE 
888:  
889:         INTEGER, INTENT(IN) :: IDI, IDJ 
890:         DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS) 
891:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_I(3), GRAD_J(3), ENERGY 
892:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_BOX(6) 
893:         LOGICAL, INTENT(IN) :: GTEST 
894:         TYPE(BOX_DERIV), INTENT(IN) :: BD 
895:         DOUBLE PRECISION, INTENT(INOUT), OPTIONAL :: GRAD_PARAMS(3) 
896:  
897:         ! RJ: particle coordinates in crystallographic space 
898:         ! RIJ: particle-particle vector in crystallographic space 
899:         ! YIJ: particle-particle vector in orthogonal space 
900:         ! MODYIJ: particle-particle distance in orthogonal space 
901:         ! DVDY: derivative of pair potential wrt MODYIJ 
902:         ! TEMP_GRAD: temporary gradient store 
903:         DOUBLE PRECISION :: RJ(3), RIJ(3), YIJ(3), MODYIJ 
904:         DOUBLE PRECISION :: DVDY, TEMP_GRAD(6) 
905:  
906:         ! PER_k: integers for looping over cell possibilities 
907:         ! CELL_RANGE: number of cells in each direction required 
908:         ! I, J: loop indices 
909:         INTEGER :: PER_X, PER_Y, PER_Z 
910:         INTEGER :: I, J 
911:  
912:         ! Loop over different coordinate possibilities. 
913:         DO PER_X = - BD%CELL_RANGE(1), BD%CELL_RANGE(1) 
914:             DO PER_Y = - BD%CELL_RANGE(2), BD%CELL_RANGE(2) 
915:                 DO PER_Z = - BD%CELL_RANGE(3), BD%CELL_RANGE(3) 
916:  
917:                     ! Skip the interaction of a particle with itself in the same 
918:                     ! cell and don't double count interactions within the unit 
919:                     ! cell. 
920:                     IF (PER_X .EQ. 0 .AND. PER_Y .EQ. 0 .AND. PER_Z .EQ. 0 & 
921:                         .AND. IDI .EQ. IDJ) CYCLE 
922:  
923:                     ! Adjust the j coordinates for the periodic image 
924:                     RJ(1) = COORDS(3*IDJ-2) + PER_X 
925:                     RJ(2) = COORDS(3*IDJ-1) + PER_Y 
926:                     RJ(3) = COORDS(3*IDJ  ) + PER_Z 
927:  
928:                     ! Find the distance in crystallographic space 
929:                     RIJ(:) = (COORDS(3*IDI-2:3*IDI) - RJ(:)) 
930:  
931:                     ! Find the distance in orthogonal space 
932:                     YIJ(:) = MATMUL(BD%ORTHOG(:,:), RIJ(:)) 
933:                     MODYIJ = SQRT(DOT_PRODUCT(YIJ(:), YIJ(:))) 
934:  
935:                     IF (MODYIJ .LT. OPP_RCUT) THEN 
936:                         ! Add energy and gradients 
937:                         ! Divide by 2 for images of the same atom, to avoid 
938:                         ! double counting 
939:                         ENERGY = ENERGY + MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * & 
940:                                           CALC_ENERGY(MODYIJ) 
941:  
942:                         IF (GTEST) THEN 
943:                             ! Divide gradient by 2 for images of the same atom, 
944:                             ! to avoid double counting 
945:                             DVDY = MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * & 
946:                                    CALC_GRAD(MODYIJ) 
947:  
948:                             ! Atom position derivatives 
949:                             DO I = 1, 3 
950:                                 TEMP_GRAD(I) = DVDY / MODYIJ * & 
951:                                              DOT_PRODUCT(YIJ(:), BD%ORTHOG(:,I)) 
952:                             END DO 
953:                             GRAD_I(:) = GRAD_I(:) + TEMP_GRAD(1:3) 
954:                             GRAD_J(:) = GRAD_J(:) - TEMP_GRAD(1:3) 
955:  
956:                             ! Lattice parameter derivatives 
957:                             TEMP_GRAD(:) = 0.D0 
958:                             DO I = 1, 6 ! Loop over the different box parameters 
959:                                 DO J = 1, 3 ! Loop over x, y, z contributions 
960:                                     TEMP_GRAD(I) = TEMP_GRAD(I) + & 
961:                                     YIJ(J)*DOT_PRODUCT(BD%DERIV(J,:,I), RIJ(:)) 
962:                                 END DO 
963:                             END DO 
964:                             TEMP_GRAD(:) = TEMP_GRAD(:) * DVDY / MODYIJ 
965:                             GRAD_BOX(:) = GRAD_BOX(:) + TEMP_GRAD(:) 
966:  
967:                             ! If necessary, calculate the parameter derivatives 
968:                             ! Divide gradient by 2 for images of the same atom, 
969:                             ! to avoid double counting 
970:                             IF(PRESENT(GRAD_PARAMS) .AND. OPP_MODE .EQ. 2)& 
971:                                 GRAD_PARAMS(:) = GRAD_PARAMS(:) + & 
972:                                     MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * & 
973:                                     CALC_GRAD_PARAMS(MODYIJ) 
974:  
975:                         END IF ! GTEST 
976:                     END IF ! End if less than cutoff 
977:                 END DO ! z loop 
978:             END DO ! y loop 
979:         END DO ! x loop 
980:  
981:     END SUBROUTINE OPP_PER_TRI 
982:  
983: !------------------------------------------------------------------------------- 
984: ! 
985: ! Calculates the contribution to the derivatives of the potential wrt the 
986: ! potential parameters OPP_K and OPP_PHI 
987: ! MODYIJ: distance between the particles for this contribution 
988: ! 
989: !------------------------------------------------------------------------------- 
990:  
991:     PURE FUNCTION CALC_GRAD_PARAMS(MODYIJ) 
992:  
993:         IMPLICIT NONE 
994:  
995:         DOUBLE PRECISION :: CALC_GRAD_PARAMS(3) 
996:         DOUBLE PRECISION, INTENT(IN) :: MODYIJ 
997:  
998:         ! S, C: sin and cos factor that appears repeatedly 
999:         ! T: working store for the output 
1000:         DOUBLE PRECISION :: S, C 
1001:         DOUBLE PRECISION :: T(3) 
1002:  
1003:         ! Initialise output 
1004:         T(:) = 0.D0 
1005:  
1006:         ! Calculate factors 
1007:         S = SIN(OPP_K * (OPP_RCUT - 1) - OPP_PHI) 
1008:         C = COS(OPP_K * (OPP_RCUT - 1) - OPP_PHI) 
1009:  
1010:         IF (BTEST(OPP_PARAMS, 0)) THEN 
1011:             ! OPP_K derivative 
1012:             ! Unmodified potential term 
1013:             T(1) = - SIN(OPP_K * (MODYIJ-1) - OPP_PHI) * (MODYIJ-1) / MODYIJ**3 
1014:             ! Term to make the potential go to zero at OPP_RCUT 
1015:             T(1) = T(1) + S * (OPP_RCUT - 1) / OPP_RCUT **3 
1016:             ! Terms to make the first derivative go to zero at OPP_RCUT 
1017:             T(1) = T(1) + (MODYIJ - OPP_RCUT) * & 
1018:                 (C * OPP_K * (OPP_RCUT - 1) + S  - & 
1019:                  S * 3.D0 * (OPP_RCUT - 1) / OPP_RCUT ) / OPP_RCUT**3 
1020:             ! Terms to make the second derivative go to zero at OPP_RCUT 
1021: !            T(1) = T(1) - 0.5D0 * (MODYIJ - OPP_RCUT)**2 * & 
1022: !                   (S * OPP_K**2 * (OPP_RCUT - 1) / OPP_RCUT**3 - & 
1023: !                    C * 2.D0 * OPP_K / OPP_RCUT**3 + & 
1024: !                    C * 6.D0 * OPP_K * (OPP_RCUT - 1) / OPP_RCUT**4 + & 
1025: !                    S * 6.D0 / OPP_RCUT**4 - & 
1026: !                    S * 12.D0 * (OPP_RCUT - 1) / OPP_RCUT**5) 
1027:         END IF 
1028:         IF (BTEST(OPP_PARAMS, 1)) THEN 
1029:             ! OPP_PHI derivative 
1030:             ! Unmodified potential term 
1031:             T(2) = SIN(OPP_K * (MODYIJ - 1) - OPP_PHI) / MODYIJ**3 
1032:             ! Term to make the potential go to zero at OPP_RCUT 
1033:             T(2) = T(2) - S / OPP_RCUT **3 
1034:             ! Terms to make the first derivative go to zero at OPP_RCUT 
1035:             T(2) = T(2) + (MODYIJ - OPP_RCUT) * & 
1036:                    ( - C * OPP_K + S * 3.D0 / OPP_RCUT ) / OPP_RCUT**3 
1037:             ! Terms to make the second derivative go to zero at OPP_RCUT 
1038: !            T(2) = T(2) - 0.5D0 * (MODYIJ - OPP_RCUT)**2 * & 
1039: !                   ( - S * OPP_K**2 / OPP_RCUT**3 - & 
1040: !                    C * 6.D0 * OPP_K / OPP_RCUT**4 + S * 12.D0 / OPP_RCUT**5)            
1041:         END IF 
1042:  
1043:         CALC_GRAD_PARAMS(:) = T(:) 
1044:  
1045:     END FUNCTION CALC_GRAD_PARAMS 
1046:  
1047: !------------------------------------------------------------------------------- 
1048: ! 
1049: ! Repel the unit cell volume from approaching 0, which repels the unit cell from 
1050: ! adopting physically impossible angle combinations. Uses a WCA potential 
1051: ! ENERGY: energy of the system 
1052: ! GRAD_ANGLES: gradients for the box angles 
1053: ! BD: information on the box and its derivatives 
1054: ! GTEST: whether to calculate gradients 
1055: ! 
1056: !------------------------------------------------------------------------------- 
1057:  
1058:     SUBROUTINE CONSTRAIN_VOLUME(ENERGY, GRAD_ANGLES, BD, GTEST) 
1059:  
1060:         IMPLICIT NONE 
1061:  
1062:         DOUBLE PRECISION, INTENT(INOUT) :: ENERGY, GRAD_ANGLES(3) 
1063:         TYPE(BOX_DERIV), INTENT(IN) :: BD 
1064:         LOGICAL, INTENT(IN) :: GTEST 
1065:  
1066:         ! Add a purely repulsive (away from zero) WCA energy term 
1067:         IF (BD%V .LT. V_SIGMA**(1.D0/6.D0)) THEN 
1068:             ENERGY = ENERGY + 4.D0 * V_EPS * & 
1069:                      ((V_SIGMA / BD%V)**(12) - (V_SIGMA / BD%V)**(6)) + V_EPS 
1070:  
1071:             IF (GTEST) THEN 
1072:                 ! Add the gradients, if necessary 
1073:                 GRAD_ANGLES(:) = GRAD_ANGLES(:) + 24.D0 * V_EPS / V_SIGMA * & 
1074:                     ((V_SIGMA / BD%V)**(7) - 2.D0 * (V_SIGMA / BD%V)**(13)) * & 
1075:                     BD%V_DERIV(:) 
1076:             END IF 
1077:         END IF 
1078:  
1079:     END SUBROUTINE CONSTRAIN_VOLUME 
1080:  
1081: !------------------------------------------------------------------------------- 
1082: ! 
1083: ! Impose limits on the allowed range of the potential parameters OPP_K by adding 
1084: ! a harmonic repulsion outside the allowed range. 
1085: ! ENERGY: energy of the system 
1086: ! GRAD_PARAMS: gradients for the potential parameters 
1087: ! GTEST: whether to calculate gradients 
1088: ! 
1089: !------------------------------------------------------------------------------- 
1090:     SUBROUTINE CONSTRAIN_PARAMETERS(ENERGY, GRAD_PARAMS, GTEST) 
1091:  
1092:         IMPLICIT NONE 
1093:  
1094:         DOUBLE PRECISION, INTENT(INOUT) :: ENERGY, GRAD_PARAMS(3) 
1095:         LOGICAL, INTENT(IN) :: GTEST 
1096:  
1097:         ! EN_CON: constraint contribution to the energy 
1098:         DOUBLE PRECISION :: EN_CON 
1099:  
1100:         ! Initialise output variable 
1101:         EN_CON = 0.D0 
1102:  
1103:         ! OPP_K energy contribution 
1104:         IF (BTEST(OPP_PARAMS, 0)) THEN 
1105:             EN_CON = EN_CON + 0.5D0 * K_REPULSION * & 
1106:                 (MERGE(OPP_K - K_UPPER, 0.D0, & 
1107:                        OPP_K .GT. K_UPPER)**2 + & 
1108:                  MERGE(OPP_K - K_LOWER, 0.D0, & 
1109:                        OPP_K .LT. K_LOWER)**2) 
1110:             IF (GTEST) THEN 
1111:             ! OPP_K gradient contribution 
1112:                 GRAD_PARAMS(1) = GRAD_PARAMS(1) + K_REPULSION * & 
1113:                     (MERGE(OPP_K - K_UPPER, 0.D0, & 
1114:                            OPP_K .GT. K_UPPER) + & 
1115:                      MERGE(OPP_K - K_LOWER, 0.D0, & 
1116:                            OPP_K .LT. K_LOWER)) 
1117:             END IF 
1118:         END IF 
1119:  
1120:         ! Add the energy contribution 
1121:         ENERGY = ENERGY + EN_CON 
1122:  
1123:     END SUBROUTINE CONSTRAIN_PARAMETERS 
1124:  
1125: !------------------------------------------------------------------------------- 
1126: ! 
1127: ! Outputs the coordinates to the file ljgauss.xyz 
1128: ! 
1129: !------------------------------------------------------------------------------- 
1130:     SUBROUTINE VIEW_OPP() 
1131:  
1132:         ! NATOMS: number of coordinates 
1133:         ! NSAVE: number of minima to write 
1134:         ! BOXx: Box dimensions 
1135:         ! PERIODIC: whether periodic boundary conditions are in use 
1136:         USE COMMONS, ONLY: NATOMS, NSAVE, PERIODIC, BOXLX, BOXLY, BOXLZ 
1137:  
1138:         ! QMIN: energies of the saved minima 
1139:         ! QMINP: coordinates of the saved minima 
1140:         ! FF: step at which the saved minima was first found 
1141:         USE QMODULE, ONLY: QMIN, QMINP, FF 
1142:  
1143:         IMPLICIT NONE 
1144:  
1145:         ! GETUNIT: function for fiding a free file unit 
1146:         ! FUNIT: output file unit 
1147:         ! NMOL: actual number of particles 
1148:         ! I, J: loop iterators 
1149:         INTEGER :: GETUNIT, FUNIT, NMOL, I, J 
1150:  
1151:         ! BOX_SIZE: convenience array for the box sizes 
1152:         ! BOX_ANGLES: convenience array for the box angles 
1153:         ! COORD: coordinate in orthogonal space 
1154:         DOUBLE PRECISION :: BOX_SIZE(3), BOX_ANGLES(3), COORD(3) 
1155:  
1156:         ! Information to transform from crystallographic to orthogonal space 
1157:         TYPE(BOX_DERIV) :: BD 
1158:  
1159:         ! Set NMOL and box parameters 
1160:         SELECT CASE (OPP_MODE) 
1161:         CASE (0) 
1162:             NMOL = NATOMS 
1163:             IF (PERIODIC) THEN 
1164:                 BOX_SIZE(1) = BOXLX 
1165:                 BOX_SIZE(2) = BOXLY 
1166:                 BOX_SIZE(3) = BOXLZ 
1167:             ELSE 
1168:                 BOX_SIZE(:) = 1.D0 
1169:             END IF 
1170:             BOX_ANGLES(:) = PI/2.D0 
1171:         CASE (1) 
1172:             NMOL = NATOMS - 2 
1173:         CASE (2) 
1174:             NMOL = NATOMS - 3 
1175:         END SELECT 
1176:  
1177:         ! Open the file 
1178:         FUNIT = GETUNIT() 
1179:         OPEN(UNIT=FUNIT, FILE='opp.xyz', STATUS='UNKNOWN') 
1180:  
1181:         ! Loop over the configurations 
1182:         DO I = 1, NSAVE 
1183:  
1184:             WRITE(FUNIT,'(I6)') NMOL 
1185:  
1186:             WRITE(FUNIT, '(A,I6,A,F20.10,A,I8,A,F20.10)') 'Energy of minimum ',& 
1187:                 I, '=', QMIN(I), ' first found at step ', FF(I), & 
1188:                 ' Energy per particle = ', QMIN(I)/NMOL 
1189:  
1190:             ! Set box parameters for varying box 
1191:             IF (OPP_MODE .NE. 0) THEN 
1192:                 BOX_SIZE(:) = QMINP(I, 3*NATOMS-2:3*NATOMS) 
1193:             END IF 
1194:             IF (OPP_MODE .EQ. 1 .OR. OPP_MODE .EQ. 2) THEN 
1195:                 BOX_ANGLES(:) = QMINP(I, 3*NATOMS-5:3*NATOMS-3) 
1196:                 BD = CALC_BOX_DERIV(BOX_ANGLES(:), BOX_SIZE(:)) 
1197:             END IF 
1198:  
1199:             DO J = 1, NMOL 
1200:  
1201:                 ! If necessary, transform the coordinate from crytallographic 
1202:                 ! space to orthogonal space 
1203:                 SELECT CASE (OPP_MODE) 
1204:                 CASE (0) 
1205:                     COORD(:) = QMINP(I, 3*J-2:3*J) 
1206:                 CASE (1, 2) 
1207:                     COORD(:) = MATMUL(BD%ORTHOG(:,:), QMINP(I, 3*J-2:3*J)) 
1208:                 END SELECT 
1209:  
1210:                 IF (PERIODIC .AND. J .EQ. 1) THEN 
1211:                     WRITE(FUNIT, '(A,3F20.10,A,3(A,F20.10))') 'O ', & 
1212:                         COORD(:), ' bbox_xyz', & 
1213:                         ' 0.0', BOX_SIZE(1), ' 0.0', BOX_SIZE(2), & 
1214:                         ' 0.0', BOX_SIZE(3) 
1215:                 ELSE 
1216:                     WRITE(FUNIT, '(A,3F20.10)') 'O ', COORD(:) 
1217:                 END IF 
1218:  
1219:             END DO ! Loop over particles 
1220:  
1221:         END DO ! Loop over the configurations 
1222:  
1223:         CLOSE(FUNIT) 
1224:  
1225:     END SUBROUTINE VIEW_OPP 
1226:  
1227: END MODULE OPP_MOD 


r32795/potential.f90 2017-06-14 20:30:13.451949038 +0100 r32794/potential.f90 2017-06-14 20:30:15.711978650 +0100
 19: SUBROUTINE POTENTIAL(X, GRAD, EREAL, GRADT, SECT) 19: SUBROUTINE POTENTIAL(X, GRAD, EREAL, GRADT, SECT)
 20: ! This subroutine acts as a wrapper for all of the potentials implemented in GMIN. 20: ! This subroutine acts as a wrapper for all of the potentials implemented in GMIN.
 21: ! The potential is chosen in the 'data' file (default is Lennard-Jones particles). 21: ! The potential is chosen in the 'data' file (default is Lennard-Jones particles).
 22: ! 22: !
 23: ! In addition, this routine also applies fields, if appropriate. 23: ! In addition, this routine also applies fields, if appropriate.
 24:    USE COMMONS 24:    USE COMMONS
 25:    USE ORBITALS_MOD 25:    USE ORBITALS_MOD
 26:    USE GENRIGID, ONLY: RIGIDINIT, NRIGIDBODY, DEGFREEDOMS, ATOMRIGIDCOORDT, & 26:    USE GENRIGID, ONLY: RIGIDINIT, NRIGIDBODY, DEGFREEDOMS, ATOMRIGIDCOORDT, &
 27:                        AACONVERGENCET, AACONVERGENCE, TRANSFORMRIGIDTOC, & 27:                        AACONVERGENCET, AACONVERGENCE, TRANSFORMRIGIDTOC, &
 28:                        TRANSFORMGRAD, TRANSFORMHESSIAN, TRANSFORMCTORIGID, & 28:                        TRANSFORMGRAD, TRANSFORMHESSIAN, TRANSFORMCTORIGID, &
 29:                        GENRIGID_DISTANCECHECK, GENRIGID_COMPRESS!, AACONVERGENCE2 29:                        GENRIGID_DISTANCECHECK, GENRIGID_COMPRESS, aaconvergence2
 30:    USE MULTIPOT, ONLY: MULTIPOT_CALL 30:    USE MULTIPOT, ONLY: MULTIPOT_CALL
 31: !   USE QMODULE, ONLY: QMIN, QMINP 31:    USE QMODULE, ONLY: QMIN, QMINP
 32:    USE PERMU, ONLY: FIN 32:    USE PERMU, ONLY: FIN
 33:    USE PORFUNCS 33:    USE PORFUNCS
 34:    USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_ENERGY_AND_GRADIENT, POT_ENE_REC_C,& 34:    USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_ENERGY_AND_GRADIENT, POT_ENE_REC_C,&
 35: &                                   AMBER12_NUM_HESS 35: &                                   AMBER12_NUM_HESS
 36:    USE MODHESS, ONLY: HESS, RBAANORMALMODET 36:    USE MODHESS, ONLY: HESS, RBAANORMALMODET
 37:    USE MODAMBER9, ONLY: ATMASS1 37:    USE MODAMBER9, ONLY: ATMASS1
 38:    USE OPEP_INTERFACE_MOD, ONLY: OPEP_ENERGY_AND_GRADIENT, OPEP_NUM_HESS 38:    USE OPEP_INTERFACE_MOD, ONLY: OPEP_ENERGY_AND_GRADIENT, OPEP_NUM_HESS
 39:    USE POLIRMOD, ONLY: POLIR 39:    USE POLIRMOD, ONLY: POLIR
 40:    USE MBPOLMOD, ONLY: MBPOL 40:    USE MBPOLMOD, ONLY: MBPOL
 41:    USE SWMOD, ONLY: SWTYPE 41:    USE SWMOD, ONLY: SWTYPE
 42:    USE MODCUDALBFGS, ONLY: CUDA_ENEGRAD_WRAPPER, CUDA_NUMERICAL_HESS 42:    USE MODCUDALBFGS, ONLY: CUDA_ENEGRAD_WRAPPER, CUDA_NUMERICAL_HESS
 43:    USE GAUSS_MOD, ONLY: GFIELD 43:    USE GAUSS_MOD, ONLY: GFIELD
 44:    USE RAD_MOD, ONLY: RAD, RADR 44:    USE RAD_MOD, ONLY: RAD, RADR
 45:    USE PREC, ONLY: INT32, REAL64 45:    USE PREC, ONLY: INT32, REAL64
 46:    USE TWIST_MOD, ONLY: TWISTT, TWIST 46:    USE TWIST_MOD, ONLY: TWISTT, TWIST
 47:    USE CONVEX_POLYHEDRA_MODULE, ONLY: CONVEX_POLYHEDRA 47:    USE CONVEX_POLYHEDRA_MODULE, ONLY: CONVEX_POLYHEDRA
 48:    USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS 48:    USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS
 49:    USE OPP_MOD, ONLY: OPP 
 50:    USE LJ_VAR_RAD_MOD, ONLY: LJ_VAR_RAD 
 51:    USE EWALD 49:    USE EWALD
 52:    IMPLICIT NONE 50:    IMPLICIT NONE
 53: ! Arguments 51: ! Arguments
 54: ! TODO: Work out intents 52: ! TODO: Work out intents
 55: ! TODO: Fix array dimensions? 53: ! TODO: Fix array dimensions?
 56:    REAL(REAL64)               :: X(*) 54:    REAL(REAL64)               :: X(*)
 57:    REAL(REAL64)               :: GRAD(*) 55:    REAL(REAL64)               :: GRAD(*)
 58:    REAL(REAL64)               :: EREAL 56:    REAL(REAL64)               :: EREAL
 59:    LOGICAL, INTENT(IN)        :: GRADT 57:    LOGICAL, INTENT(IN)        :: GRADT
 60:    LOGICAL, INTENT(IN)        :: SECT 58:    LOGICAL, INTENT(IN)        :: SECT
 61: ! Variables 59: ! Variables
 62: ! hk286 > generalised rigid body additions 60: ! hk286 > generalised rigid body additions
 63:    REAL(REAL64)               :: XCOORDS(3*NATOMS), GRADATOMS(3*NATOMS), & 61:    REAL(REAL64)               :: XCOORDS(3*NATOMS), GRADATOMS(3*NATOMS), &
 64:                                  XRIGIDCOORDS(DEGFREEDOMS), XRIGIDGRAD(DEGFREEDOMS)!, XFRAC(3*NATOMS) 62:                                  XRIGIDCOORDS(DEGFREEDOMS), XRIGIDGRAD(DEGFREEDOMS), XFRAC(3*NATOMS)
 65:    REAL(REAL64), ALLOCATABLE  :: TEMPHESS(:) 63:    REAL(REAL64), ALLOCATABLE  :: TEMPHESS(:)
 66:    REAL(REAL64)               :: GRAD1(3*NATOMS) 64:    REAL(REAL64)               :: GRAD1(3*NATOMS)
 67:    INTEGER(INT32)             :: I, J, K 65:    INTEGER(INT32)             :: I, J, K
 68:    REAL(REAL64)               :: XRIGIDHESS(DEGFREEDOMS, DEGFREEDOMS) 66:    REAL(REAL64)               :: XRIGIDHESS(DEGFREEDOMS, DEGFREEDOMS)
 69:    LOGICAL                    :: FTEST, EVAP, COMPON, YESNO, GUIDET, EVAPREJECT 67:    LOGICAL                    :: FTEST, EVAP, COMPON, YESNO, GUIDET, EVAPREJECT
 70:    INTEGER(INT32)             :: J1, J2, J3, CSMIT 68:    INTEGER(INT32)             :: J1, J2, J3, CSMIT
 71:    CHARACTER                  :: FNAME*80, DUMM*4 69:    CHARACTER                  :: FNAME*80, DUMM*4
 72:    REAL(REAL64)               :: DUMMY2, GEMAX, RMAT(3, 3), XTEMP(3*NATOMS), AVVAL, & 70:    REAL(REAL64)               :: DUMMY2, GEMAX, RMAT(3, 3), XTEMP(3*NATOMS), AVVAL, &
 73:                                  DUMMY(3*NATOMS), DIST2, QE, QX, AA(3), SAVECSMNORM, & 71:                                  DUMMY(3*NATOMS), DIST2, QE, QX, AA(3), SAVECSMNORM, &
 74:                                  CSMRMS, CSMGRAD(3), SAVECSMPMAT(3, 3), & 72:                                  CSMRMS, CSMGRAD(3), SAVECSMPMAT(3, 3), &
 75:                                  SAVECSMIMAGES(3*NATOMS*CSMGPINDEX)!, BOXPARAMS(6) 73:                                  SAVECSMIMAGES(3*NATOMS*CSMGPINDEX), BOXPARAMS(6)
 76:    CHARACTER(LEN=3)           :: CSMGPSAVE 74:    CHARACTER(LEN=3)           :: CSMGPSAVE
 77:    INTEGER(INT32)             :: CSMGPINDEXSAVE 75:    INTEGER(INT32)             :: CSMGPINDEXSAVE
 78:    REAL(REAL64)               :: PTGPSAVE(3, 3, 2*CSMGPINDEX), CSMNORMSAVE, ENERGY, VNEW(3*NATOMS) 76:    REAL(REAL64)               :: PTGPSAVE(3, 3, 2*CSMGPINDEX), CSMNORMSAVE, ENERGY, VNEW(3*NATOMS)
 79: !   INTEGER(INT32)             :: BRUN 77:    INTEGER(INT32)             :: BRUN
 80:    INTEGER(INT32)             :: NQTOT 78:    INTEGER(INT32)             :: NQTOT
 81:    LOGICAL                    :: SOCOUPLE, GUIDECHANGET, CSMDOGUIDET, COMTEST 79:    LOGICAL                    :: SOCOUPLE, GUIDECHANGET, CSMDOGUIDET, COMTEST
 82: ! khs26 > AMBER12 energy decomposition, type defined in AMBER12 interface 80: ! khs26 > AMBER12 energy decomposition, type defined in AMBER12 interface
 83:    TYPE(POT_ENE_REC_C)        :: AMBER12_ENERGY_DECOMP 81:    TYPE(POT_ENE_REC_C)        :: AMBER12_ENERGY_DECOMP
 84: !Used only in commented sections 82: !Used only in commented sections
 85:    REAL(REAL64)               :: GRADDUM(3*NATOMS), EPLUS, EMINUS, DIFF!, TERMLJ, TERMMF, & 83:    REAL(REAL64)               :: XG, YG, ZG, GRADLJ(3*NATOMS), EREALLJ, GRADMF(3*NATOMS), EREALMF, TERMLJ, TERMMF, &
 86:                                  !, XG, YG, ZG, RMS1, RMS2, GRADLJ(3*NATOMS), GRADMF(3*NATOMS), & 84:                                  GRADDUM(3*NATOMS), GRADDUM1(3*NATOMS), GRADDUM2(3*NATOMS), EPLUS, EMINUS, DIFF, &
 87: !                                GRADDUM1(3*NATOMS), GRADDUM2(3*NATOMS), EREALLJ, EREALMF, EDUM1, EDUM2 85:                                  EDUM1, EDUM2, RMS1, RMS2
 88: !DS656> GRADIENT/HESSIAN TESTING 86: !DS656> GRADIENT/HESSIAN TESTING
 89:    DOUBLE PRECISION :: GPLUS(3*NATOMS), GMINUS(3*NATOMS) 87:    DOUBLE PRECISION :: GPLUS(3*NATOMS), GMINUS(3*NATOMS)
 90:  88: 
 91: ! TODO: Delete common blocks 89: ! TODO: Delete common blocks
 92:    COMMON /CO/ COMPON 90:    COMMON /CO/ COMPON
 93:    COMMON /FAIL/ FTEST 91:    COMMON /FAIL/ FTEST
 94:    COMMON /EV/ EVAP, EVAPREJECT 92:    COMMON /EV/ EVAP, EVAPREJECT
 95:    COMMON /GD/ GUIDECHANGET, GUIDET, CSMDOGUIDET 93:    COMMON /GD/ GUIDECHANGET, GUIDET, CSMDOGUIDET
 96:    COMMON /CSMAVVAL/ AVVAL, CSMRMS, CSMIT 94:    COMMON /CSMAVVAL/ AVVAL, CSMRMS, CSMIT
 97:    COMMON /TOT/ NQTOT 95:    COMMON /TOT/ NQTOT
350: !         X(1:DEGFREEDOMS) = XRIGIDCOORDS(1:DEGFREEDOMS)348: !         X(1:DEGFREEDOMS) = XRIGIDCOORDS(1:DEGFREEDOMS)
351: !      ENDIF349: !      ENDIF
352: !      DO J1=1, degfreedoms350: !      DO J1=1, degfreedoms
353: !         !print *, 'index: ', j1351: !         !print *, 'index: ', j1
354: !         X(J1)=X(J1)+DIFF352: !         X(J1)=X(J1)+DIFF
355: !         CALL benzgenrigidewald(X, GRADDUM, EPLUS,.FALSE.)353: !         CALL benzgenrigidewald(X, GRADDUM, EPLUS,.FALSE.)
356: !        !print *, 'subtracting now'354: !        !print *, 'subtracting now'
357: !         X(J1)=X(J1)-2.0D0*DIFF355: !         X(J1)=X(J1)-2.0D0*DIFF
358: !         CALL benzgenrigidewald(X, GRADDUM, EMINUS,.FALSE.)356: !         CALL benzgenrigidewald(X, GRADDUM, EMINUS,.FALSE.)
359: !         X(J1)=X(J1)+DIFF357: !         X(J1)=X(J1)+DIFF
360: !         print *, j1, 'diff', 100.0d0*abs(grad(j1) - ((eplus-eminus)/(2.0d0*diff)))/abs(grad(j1)), &358: !         print *, j1, 'diff', 100.0d0*abs(grad(j1) - ((eplus-eminus)/(2.0d0*diff)))/abs(grad(j1)), (grad(j1) - (eplus-eminus)/(2.0d0*diff))
361: !                  (grad(j1) - (eplus-eminus)/(2.0d0*diff)) 
362: !         IF ((ABS(GRAD(J1)).NE.0.0D0).AND.((100.0D0*(ABS(GRAD(J1)-((EPLUS-EMINUS)/(2.0D0*DIFF)))/ABS(GRAD(J1)))).GT.1.0D0)) THEN359: !         IF ((ABS(GRAD(J1)).NE.0.0D0).AND.((100.0D0*(ABS(GRAD(J1)-((EPLUS-EMINUS)/(2.0D0*DIFF)))/ABS(GRAD(J1)))).GT.1.0D0)) THEN
363: !             WRITE(*,'(I5, 2F20.10)') J1, GRAD(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)360: !             WRITE(*,'(I5, 2F20.10)') J1, GRAD(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)
364: !         else361: !         else
365: !            print *, 'fine', j1, grad(j1), (eplus-eminus)/(2.0d0*diff)362: !            print *, 'fine', j1, grad(j1), (eplus-eminus)/(2.0d0*diff)
366: !         END IF363: !         END IF
367: !      ENDDO364: !      ENDDO
368: !      stop365: !      stop
369: 366: 
370: ! dj337: Ewald sums367: ! dj337: Ewald sums
371:    ELSE IF (EWALDT) THEN368:    ELSE IF (EWALDT) THEN
1033: ! ds656> Multicomponent LJ with no cutoff (akin to BLJ_CLUST).1030: ! ds656> Multicomponent LJ with no cutoff (akin to BLJ_CLUST).
1034: !        Should give same results as GLJ, but different code.1031: !        Should give same results as GLJ, but different code.
1035:       CALL RAD(X, GRAD, EREAL, GRADT)1032:       CALL RAD(X, GRAD, EREAL, GRADT)
1036:       CALL MLJ(X, GRAD, EREAL, GRADT, SECT, .FALSE.)1033:       CALL MLJ(X, GRAD, EREAL, GRADT, SECT, .FALSE.)
1037: 1034: 
1038:    ELSE IF (GLJT) THEN1035:    ELSE IF (GLJT) THEN
1039: ! generalised LJ with no cutoff1036: ! generalised LJ with no cutoff
1040:       CALL RAD(X, GRAD, EREAL, GRADT)1037:       CALL RAD(X, GRAD, EREAL, GRADT)
1041:       CALL GLJ(X, GRAD, EREAL, GRADT)1038:       CALL GLJ(X, GRAD, EREAL, GRADT)
1042: 1039: 
1043:    ELSE IF (LJ_VAR_RADT) THEN 
1044: ! jwrm2> Lennard-Jones with each particle having a different radius. Works with 
1045: !        RIGIDINIT 
1046:       IF (RIGIDINIT .AND. .NOT. ATOMRIGIDCOORDT) THEN 
1047:          XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS) 
1048:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS) 
1049:       END IF 
1050:       CALL LJ_VAR_RAD(X, GRAD, EREAL, GRADT) 
1051:       IF (RIGIDINIT .AND. .NOT. ATOMRIGIDCOORDT) THEN 
1052:          X(DEGFREEDOMS+1:3*NATOMS)=0.0D0 
1053:          X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS) 
1054:          IF (GRADT) THEN 
1055:             GRADATOMS(1:3*NATOMS) = GRAD(1:3*NATOMS) ! Needed for AACONVERGENCE 
1056:             CALL TRANSFORMGRAD(GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD) 
1057:             GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0 
1058:             GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS) 
1059:          END IF 
1060:      END IF 
1061:  
1062:    ELSE IF (BINARY) THEN1040:    ELSE IF (BINARY) THEN
1063:       IF (SHIFTCUT) THEN1041:       IF (SHIFTCUT) THEN
1064:          if (rigidinit.and.(.not.atomrigidcoordt)) then1042:          if (rigidinit.and.(.not.atomrigidcoordt)) then
1065:             xrigidcoords(1:degfreedoms) = x(1:degfreedoms)1043:             xrigidcoords(1:degfreedoms) = x(1:degfreedoms)
1066:             call transformrigidtoc(1, nrigidbody, x, xrigidcoords)1044:             call transformrigidtoc(1, nrigidbody, x, xrigidcoords)
1067:          endif1045:          endif
1068:          call ljpshift(x, grad, ereal, gradt, sect)1046:          call ljpshift(x, grad, ereal, gradt, sect)
1069:          if (rigidinit.and.(.not.atomrigidcoordt)) then1047:          if (rigidinit.and.(.not.atomrigidcoordt)) then
1070:             x(degfreedoms+1:3*natoms)=0.0d01048:             x(degfreedoms+1:3*natoms)=0.0d0
1071:             x(1:degfreedoms)=xrigidcoords(1:degfreedoms)1049:             x(1:degfreedoms)=xrigidcoords(1:degfreedoms)
1087:       CALL RADR(X, GRAD, EREAL, GRADT)1065:       CALL RADR(X, GRAD, EREAL, GRADT)
1088:       CALL LB2(X, GRAD, EREAL, GRADT)1066:       CALL LB2(X, GRAD, EREAL, GRADT)
1089: 1067: 
1090:    ELSE IF (LJCOULT) THEN1068:    ELSE IF (LJCOULT) THEN
1091:       CALL RADR(X, GRAD, EREAL, GRADT)1069:       CALL RADR(X, GRAD, EREAL, GRADT)
1092:       CALL LJCOUL(X, GRAD, EREAL, GRADT)1070:       CALL LJCOUL(X, GRAD, EREAL, GRADT)
1093: 1071: 
1094:   ELSE IF (LJ_GAUSST) THEN1072:   ELSE IF (LJ_GAUSST) THEN
1095:       CALL LJ_GAUSS(X, GRAD, EREAL, GRADT)1073:       CALL LJ_GAUSS(X, GRAD, EREAL, GRADT)
1096: 1074: 
1097:   ELSE IF (OPPT) THEN 
1098:       CALL OPP(X, GRAD, EREAL, GRADT) 
1099:  
1100:    ELSE IF (STOCKT) THEN1075:    ELSE IF (STOCKT) THEN
1101:       CALL RADR(X, GRAD, EREAL, GRADT)1076:       CALL RADR(X, GRAD, EREAL, GRADT)
1102:       CALL STOCK(X, GRAD, EREAL, GRADT)1077:       CALL STOCK(X, GRAD, EREAL, GRADT)
1103: 1078: 
1104: ! RBAA potentials1079: ! RBAA potentials
1105: 1080: 
1106:    ELSE IF (CAPBINT) THEN1081:    ELSE IF (CAPBINT) THEN
1107:       CALL CAPBIN (X, GRAD, EREAL, GRADT)1082:       CALL CAPBIN (X, GRAD, EREAL, GRADT)
1108: 1083: 
1109:    ELSE IF (CHIROT) THEN1084:    ELSE IF (CHIROT) THEN


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0