hdiff output

r33077/commons.f90 2017-07-26 19:30:27.876719238 +0100 r33076/commons.f90 2017-07-26 19:30:30.008747549 +0100
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, READMASST, SPECMASST, NEWTSALLIST, &117:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, LJ_VAR_RADT, 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, RIGIDMBPOLT, &119:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, RIGIDMBPOLT, &
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, EWALDT, WATERMETHANET, MLPVB3NNT, CLATHRATET124:      &        SQNMT, SQNM_DEBUGT, SQNM_BIOT, BENZRIGIDEWALDT, ORTHO, EWALDT, WATERMETHANET, MLPVB3NNT, CLATHRATET
125: !125: !
126:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:)126:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:)
127:       DOUBLE PRECISION, ALLOCATABLE :: ATMASS(:)127:       DOUBLE PRECISION, ALLOCATABLE :: ATMASS(:)


r33077/GMIN.tex 2017-07-26 19:30:27.608715678 +0100 r33076/GMIN.tex 2017-07-26 19:30:29.752744150 +0100
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, 2 or 3. \textit{mode} 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. \textit{mode} 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)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,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.1135: and if the 4 bit is set, $\sigma^2$ will be optimised.
1136: 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: 
1137: \item {\it LOCALSAMPLE abthresh acthresh\/}: Keyword currently under construction! For three groups of atoms defined in movableatoms1141: \item {\it LOCALSAMPLE abthresh acthresh\/}: Keyword currently under construction! For three groups of atoms defined in movableatoms
1138: (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}. 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}. 
1139: If this condition is broken AFTER the quench, it is automatically rejected.1143: If this condition is broken AFTER the quench, it is automatically rejected.
1140: 1144: 
1141: \item {\it LPERMDIST n thresh cut orbittol\/}: provides the preferred local1145: \item {\it LPERMDIST n thresh cut orbittol\/}: provides the preferred local
1142: permutational alignment procedure.1146: permutational alignment procedure.
1143: {\it n\/} is the maximum number of extra atoms to include in running1147: {\it n\/} is the maximum number of extra atoms to include in running
1144: {\bf minperm\/} for each permutable group (default 4), rather than1148: {\bf minperm\/} for each permutable group (default 4), rather than
1145: treating the whole system at the same time, as for {\it PERMDIST}.1149: treating the whole system at the same time, as for {\it PERMDIST}.
1146: The extra atoms that define the neighbourhood of the permutable group are1150: The extra atoms that define the neighbourhood of the permutable group are
1845: {\it x1 y1 z1} \\1849: {\it x1 y1 z1} \\
1846: {\it x2 y2 z2} \\1850: {\it x2 y2 z2} \\
1847: {\it \ldots} \\1851: {\it \ldots} \\
1848: {\it \ldots} \\1852: {\it \ldots} \\
1849: {\it xn yn zn} \\1853: {\it xn yn zn} \\
1850: 1854: 
1851: The file {\it rbodyconfig\/} defines the rigid bodies, with the following format: \\1855: The file {\it rbodyconfig\/} defines the rigid bodies, with the following format: \\
1852: {\it GROUP NoAtomsInGroup} \\1856: {\it GROUP NoAtomsInGroup} \\
1853: {\it Atom 1} \\1857: {\it Atom 1} \\
1854: {\it Atom 2} \\1858: {\it Atom 2} \\
1855: {\it \ldots} \\1859: {\it ldots} \\
1856: {\it \ldots} \\1860: {\it ldots} \\
1857: {\it Atom N} \\1861: {\it Atom N} \\
1858: The keywords {\it RELAXFINALQUENCH\/} and {\it NRELAXRIGID\/} can be used in conjuction with {\it RIGIDINIT\/}.1862: The keywords {\it RELAXFINALQUENCH\/} and {\it NRELAXRIGID\/} can be used in conjuction with {\it RIGIDINIT\/}.
1859: 1863: 
1860: \item {\it RINGROTSCALE factor\/}: when applying cartesian moves with CHARMM, amino acid rings are moved as rigid units. Setting {\it factor} (default 0.0) between 0.0 and 1.0 will apply a random rotation to these rings during step taking. The suggested value is 0.1 to prevent the regular formation of high energy structures. 1864: \item {\it RINGROTSCALE factor\/}: when applying cartesian moves with CHARMM, amino acid rings are moved as rigid units. Setting {\it factor} (default 0.0) between 0.0 and 1.0 will apply a random rotation to these rings during step taking. The suggested value is 0.1 to prevent the regular formation of high energy structures. 
1861: 1865: 
1862: \item {\it RKMIN\/}: specifies a Runga-Kutta minimisation scheme. 1866: \item {\it RKMIN\/}: specifies a Runga-Kutta minimisation scheme. 
1863: Very inefficient.1867: Very inefficient.
1864: 1868: 
1865: \item{\it RMS rmslimit rmstol rmssave lca best}: used with {\it CHARMM} keyword to1869: \item{\it RMS rmslimit rmstol rmssave lca best}: used with {\it CHARMM} keyword to
1866: specify that the RMSD compared to a reference geometry is calculated. The reference geometry must 1870: specify that the RMSD compared to a reference geometry is calculated. The reference geometry must 


r33077/keywords.f 2017-07-26 19:30:28.144722797 +0100 r33076/keywords.f 2017-07-26 19:30:30.296751374 +0100
 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, 46:       USE OPP_MOD, ONLY: OPP_MODE, OPP_RCUT, OPP_K, OPP_PHI, OPP_PARAMS,
 47:      &                   OPP_INITIALISE 47:      &                   OPP_INITIALISE
  48:       USE LJ_VAR_RAD_MOD, ONLY: LJ_VAR_RAD_INITIALISE
 48:       USE MBPOLMOD, ONLY: MBPOLINIT 49:       USE MBPOLMOD, ONLY: MBPOLINIT
 49:       USE SWMOD, ONLY: SWINIT, MWINIT 50:       USE SWMOD, ONLY: SWINIT, MWINIT
 50:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_GET_COORDS 51:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_GET_COORDS
 51: !     &                                 AMBER12_ATOMSS, AMBER12_SETUP, 52: !     &                                 AMBER12_ATOMSS, AMBER12_SETUP,
 52: !     &                                 AMBER12_RESIDUES, 53: !     &                                 AMBER12_RESIDUES,
 53: !     &                                 POPULATE_ATOM_DATA 54: !     &                                 POPULATE_ATOM_DATA
 54:       USE CHIRALITY, ONLY : CIS_TRANS_TOL 55:       USE CHIRALITY, ONLY : CIS_TRANS_TOL
 55:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR 56:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR
 56:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS, 57:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS,
 57:      &     PARSE_MLJ_PARAMS 58:      &     PARSE_MLJ_PARAMS
1130: !cv320> rigid body watermethane1131: !cv320> rigid body watermethane
1131:       WATERMETHANET=.FALSE.1132:       WATERMETHANET=.FALSE.
1132:       CLATHRATET= .FALSE.1133:       CLATHRATET= .FALSE.
1133: ! ds656> Parallelised generalised basin-hopping1134: ! ds656> Parallelised generalised basin-hopping
1134:       GBHT=.FALSE.1135:       GBHT=.FALSE.
1135: 1136: 
1136: ! General mixed LJ systems1137: ! General mixed LJ systems
1137:       GLJT=.FALSE.1138:       GLJT=.FALSE.
1138: ! ds656> Multicomponent LJ system (different implementation to GLJ!)1139: ! ds656> Multicomponent LJ system (different implementation to GLJ!)
1139:       MLJT=.FALSE.1140:       MLJT=.FALSE.
 1141: ! jwrm2> LJ with each atom having a different radius. Works with RIGIDINIT.
 1142:       LJ_VAR_RADT = .FALSE.
1140: 1143: 
1141: ! khs26> Free energy basin-hopping stuff1144: ! khs26> Free energy basin-hopping stuff
1142:       FEBHT = .FALSE.1145:       FEBHT = .FALSE.
1143:       FETEMP = 0.0D01146:       FETEMP = 0.0D0
1144: ! khs26> This requires a minimum separation between the zero and non-zero1147: ! khs26> This requires a minimum separation between the zero and non-zero
1145: ! eigenvalues for runs using free energy basin-hopping. This is based on1148: ! eigenvalues for runs using free energy basin-hopping. This is based on
1146: ! the minimum found in quench zero.1149: ! the minimum found in quench zero.
1147:       MIN_ZERO_SEP = 0.0D01150:       MIN_ZERO_SEP = 0.0D0
1148:       MAX_ATTEMPTS = 51151:       MAX_ATTEMPTS = 5
1149: ! Use sparse hessian methods1152: ! Use sparse hessian methods
4843: 4846: 
4844: !4847: !
4845: ! NOT DOCUMENTED4848: ! NOT DOCUMENTED
4846: !4849: !
4847:       ELSE IF (WORD.EQ.'LJMF') THEN4850:       ELSE IF (WORD.EQ.'LJMF') THEN
4848:          LJMFT=.TRUE.4851:          LJMFT=.TRUE.
4849: !        CALL LJPARAMMF4852: !        CALL LJPARAMMF
4850:          WRITE(MYUNIT,'(A)') 'LJMF not currently maintained'4853:          WRITE(MYUNIT,'(A)') 'LJMF not currently maintained'
4851:          STOP4854:          STOP
4852: 4855: 
 4856: ! jwrm2> LJ with each atom having a different radius. Works with RIGIDINIT.
 4857:       ELSE IF (WORD .EQ. 'LJVARRAD') THEN
 4858:          LJ_VAR_RADT = .TRUE.
 4859:          CALL LJ_VAR_RAD_INITIALISE()
 4860: 
4853: ! start = an N-oligomer is constructed by relocating NREPEAT units and placing them4861: ! start = an N-oligomer is constructed by relocating NREPEAT units and placing them
4854: ! at random distance R with Rmin <= R <= Rmax and angle alpha in the xy-plane4862: ! at random distance R with Rmin <= R <= Rmax and angle alpha in the xy-plane
4855: ! transxy: rigid body translation only in the xy-plane4863: ! transxy: rigid body translation only in the xy-plane
4856: ! rotz:  rigid body rotation only around the z-axis4864: ! rotz:  rigid body rotation only around the z-axis
4857: ! dihesc: only perturbation to the side chains4865: ! dihesc: only perturbation to the side chains
4858: !4866: !
4859:       ELSE IF (WORD.EQ.'MAKEOLIGO') THEN4867:       ELSE IF (WORD.EQ.'MAKEOLIGO') THEN
4860:          MAKEOLIGOT=.TRUE.4868:          MAKEOLIGOT=.TRUE.
4861:          CALL READA(UNSTRING)4869:          CALL READA(UNSTRING)
4862:          IF (TRIM(ADJUSTL(UNSTRING)).EQ.'START') THEN4870:          IF (TRIM(ADJUSTL(UNSTRING)).EQ.'START') THEN


r33077/lj_gauss.f90 2017-07-26 19:30:28.404726249 +0100 r33076/lj_gauss.f90 2017-07-26 19:30:30.564754933 +0100
 33: PUBLIC :: LJ_GAUSS, LJ_GAUSS_INITIALISE, VIEW_LJ_GAUSS, LJ_GAUSS_TAKESTEP 33: PUBLIC :: LJ_GAUSS, LJ_GAUSS_INITIALISE, VIEW_LJ_GAUSS, LJ_GAUSS_TAKESTEP
 34: ! Private parameters 34: ! Private parameters
 35: PRIVATE :: RCUTINV6, RCUTINV12, EXP_RCUT, PREF_RCUT, PI, SCL_FCT 35: PRIVATE :: RCUTINV6, RCUTINV12, EXP_RCUT, PREF_RCUT, PI, SCL_FCT
 36: PRIVATE :: IMAGE_CUTOFF, MAX_LENGTH_STEP, MAX_ANGLE_STEP, MAX_EPS_STEP 36: PRIVATE :: IMAGE_CUTOFF, MAX_LENGTH_STEP, MAX_ANGLE_STEP, MAX_EPS_STEP
 37: PRIVATE :: MAX_R0_STEP, MAX_SIGMASQ_STEP, V_EPS, V_SIGMA, EPS_LOWER, EPS_UPPER 37: PRIVATE :: MAX_R0_STEP, MAX_SIGMASQ_STEP, V_EPS, V_SIGMA, EPS_LOWER, EPS_UPPER
 38: PRIVATE :: EPS_REPULSION, R0_LOWER, R0_UPPER, R0_REPULSION, SIGMASQ_LOWER 38: PRIVATE :: EPS_REPULSION, R0_LOWER, R0_UPPER, R0_REPULSION, SIGMASQ_LOWER
 39: PRIVATE :: SIGMASQ_UPPER, SIGMASQ_REPULSION 39: PRIVATE :: SIGMASQ_UPPER, SIGMASQ_REPULSION
 40: ! Private functions 40: ! Private functions
 41: PRIVATE :: PERIODIC_RESET, LJ_GAUSS_PER, CALC_ENERGY, CALC_ENERGY_LJ 41: PRIVATE :: PERIODIC_RESET, LJ_GAUSS_PER, CALC_ENERGY, CALC_ENERGY_LJ
 42: PRIVATE :: CALC_ENERGY_GAUSS, CALC_GRAD, CALC_FCTS, CONSTRAIN_PARAMETERS 42: PRIVATE :: CALC_ENERGY_GAUSS, CALC_GRAD, CALC_FCTS, CONSTRAIN_PARAMETERS
 43: PRIVATE :: CHECK_ANGLES, REJECT, CONSTRAIN_VOLUME 43: PRIVATE :: BOX_DERIV, CALC_BOX_DERIV, CHECK_ANGLES, REJECT, CONSTRAIN_VOLUME
 44:  44: 
 45: ! EPS: strength of the Gaussian potential. The LJ strength is 1. 45: ! EPS: strength of the Gaussian potential. The LJ strength is 1.
 46: ! R0: Equilibrium distance of the minimum in the Gaussian potential. The LJ 46: ! R0: Equilibrium distance of the minimum in the Gaussian potential. The LJ
 47: !     rmin is 1. 47: !     rmin is 1.
 48: ! SIGMASQ: variation (width) of the Gaussian potential 48: ! SIGMASQ: variation (width) of the Gaussian potential
 49: ! RCUTSQ: cutoff distance at which the energy, gradients and Hessian of 49: ! RCUTSQ: cutoff distance at which the energy, gradients and Hessian of
 50: !         both parts of the potential go smoothly to zero. 50: !         both parts of the potential go smoothly to zero.
 51: DOUBLE PRECISION :: LJ_GAUSS_EPS, LJ_GAUSS_R0, LJ_GAUSS_SIGMASQ, LJ_GAUSS_RCUT 51: DOUBLE PRECISION :: LJ_GAUSS_EPS, LJ_GAUSS_R0, LJ_GAUSS_SIGMASQ, LJ_GAUSS_RCUT
 52:  52: 
 53: ! LJ_GAUSS_MODE: 0 standard minimisation 53: ! LJ_GAUSS_MODE: 0 standard minimisation
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.D9
115: 115: 
 116: ! BOX_DERIV: stores information for calculating lattice derivatives
 117: TYPE :: BOX_DERIV
 118:     ! ORTHOG: the orthogonal tansformation matrix, for changing from
 119:     !         crystallographic space to orthogonal space first index is row,
 120:     !         second is column
 121:     DOUBLE PRECISION :: ORTHOG(3,3)
 122: 
 123:     ! DERIV: derivatives of the orthogonal matrix wrt lattice parameters,
 124:     !        first index indicates the row, second index the column and
 125:     !        third index the parameter, with 1-3 being the angles and 4-6 being
 126:     !        the lengths
 127:     DOUBLE PRECISION :: DERIV(3, 3, 6)
 128: 
 129:     ! RECIP: lengths of the reciprocal lattice vectors
 130:     DOUBLE PRECISION :: RECIP(3)
 131: 
 132:     ! CELL_RANGE: number of cells in each direction that must be evaluated
 133:     INTEGER :: CELL_RANGE(3)
 134: 
 135:     ! V: dimensionless volume of the unit cell, ie volume if all lengths were 1
 136:     DOUBLE PRECISION :: V
 137: 
 138:     ! V_DERIV: derivatievs of V wrt the lattice angles
 139:     DOUBLE PRECISION :: V_DERIV(3)
 140: END TYPE BOX_DERIV
 141: 
116: CONTAINS142: CONTAINS
117: 143: 
118: !-------------------------------------------------------------------------------144: !-------------------------------------------------------------------------------
119: !145: !
120: ! Main routine for calculating the energy and gradients.146: ! Main routine for calculating the energy and gradients.
121: ! X: coordinates array147: ! X: coordinates array
122: ! GRAD: gradients array148: ! GRAD: gradients array
123: ! ENERGY: calculated energy149: ! ENERGY: calculated energy
124: ! GTEST: whether gradients should be calculated150: ! GTEST: whether gradients should be calculated
125: !151: !
126: !-------------------------------------------------------------------------------152: !-------------------------------------------------------------------------------
127:     SUBROUTINE LJ_GAUSS(X, GRAD, ENERGY, GTEST)153:     SUBROUTINE LJ_GAUSS(X, GRAD, ENERGY, GTEST)
128: 154: 
129:         ! MYUNIT: file unit for main output file155:         ! MYUNIT: file unit for main output file
130:         ! NATOMS: number of particles156:         ! NATOMS: number of particles
131:         ! PERIODIC: whether to use periodic boundary conditions157:         ! PERIODIC: whether to use periodic boundary conditions
132:         USE COMMONS, ONLY: MYUNIT, NATOMS, PERIODIC158:         USE COMMONS, ONLY: MYUNIT, NATOMS, PERIODIC
133: 159: 
134:         ! BOX_DERIV: type for storing box derivative information 
135:         ! CALC_BOX_DERIV: calculates box derivative information 
136:         USE TRICLINIC_MOD, ONLY: BOX_DERIV, CALC_BOX_DERIV 
137:  
138:         IMPLICIT NONE160:         IMPLICIT NONE
139: 161: 
140:         DOUBLE PRECISION, INTENT(INOUT) :: X(3*NATOMS)162:         DOUBLE PRECISION, INTENT(INOUT) :: X(3*NATOMS)
141:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOMS), ENERGY163:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOMS), ENERGY
142:         LOGICAL, INTENT(IN) :: GTEST164:         LOGICAL, INTENT(IN) :: GTEST
143: 165: 
144:         ! I, J: indices for looping over particles166:         ! I, J: indices for looping over particles
145:         ! NMOL: actual number of particles167:         ! NMOL: actual number of particles
146:         INTEGER :: I, J, NMOL168:         INTEGER :: I, J, NMOL
147: 169: 
230:                     END DO ! Inner loop over particles252:                     END DO ! Inner loop over particles
231:                 END DO ! Outer loop over particles253:                 END DO ! Outer loop over particles
232: 254: 
233:             CASE(2)255:             CASE(2)
234:                 ! Triclinic varying lattice256:                 ! Triclinic varying lattice
235:                 ! Reset all particles to within the box257:                 ! Reset all particles to within the box
236:                 CALL PERIODIC_RESET(X(:))258:                 CALL PERIODIC_RESET(X(:))
237: 259: 
238:                 ! Calculate box derivative factors260:                 ! Calculate box derivative factors
239:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &261:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &
240:                                     X(3*NATOMS-2:3*NATOMS  ), LJ_GAUSS_RCUT)262:                                     X(3*NATOMS-2:3*NATOMS  ))
241: 263: 
242: !                WRITE(MYUNIT, *) 'Box lengths: ', X(3*NATOMS-2:3*NATOMS)264: !                WRITE(MYUNIT, *) 'Box lengths: ', X(3*NATOMS-2:3*NATOMS)
243: !                WRITE(MYUNIT, *) 'Box angles: ', X(3*NATOMS-5:3*NATOMS-3)265: !                WRITE(MYUNIT, *) 'Box angles: ', X(3*NATOMS-5:3*NATOMS-3)
244: !                WRITE(MYUNIT, *) 'CELL_RANGE(:) ', BD%CELL_RANGE(:)266: !                WRITE(MYUNIT, *) 'CELL_RANGE(:) ', BD%CELL_RANGE(:)
245: !                WRITE(MYUNIT, *) 'V = ', BD%V267: !                WRITE(MYUNIT, *) 'V = ', BD%V
246: 268: 
247:                 ! Check whether the unit cell angles are physically possible.269:                 ! Check whether the unit cell angles are physically possible.
248:                 ! Reject if not.270:                 ! Reject if not.
249:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN271:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN
250:                     CALL REJECT(ENERGY, GRAD)272:                     CALL REJECT(ENERGY, GRAD)
279:                 CALL CONSTRAIN_VOLUME(ENERGY, GRAD(3*NATOMS-5:3*NATOMS-3), BD, &301:                 CALL CONSTRAIN_VOLUME(ENERGY, GRAD(3*NATOMS-5:3*NATOMS-3), BD, &
280:                                       GTEST)302:                                       GTEST)
281: 303: 
282:             CASE(3)304:             CASE(3)
283:                 ! Triclinic varying lattice, with varying parameters305:                 ! Triclinic varying lattice, with varying parameters
284:                 ! Reset all particles to within the box306:                 ! Reset all particles to within the box
285:                 CALL PERIODIC_RESET(X(:))307:                 CALL PERIODIC_RESET(X(:))
286: 308: 
287:                 ! Calculate box derivative factors309:                 ! Calculate box derivative factors
288:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &310:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &
289:                                     X(3*NATOMS-2:3*NATOMS  ), LJ_GAUSS_RCUT)311:                                     X(3*NATOMS-2:3*NATOMS  ))
290: 312: 
291:                 ! Check whether the unit cell angles are physically possible.313:                 ! Check whether the unit cell angles are physically possible.
292:                 ! If we have a disallowed unit cell, we set the energy to very314:                 ! If we have a disallowed unit cell, we set the energy to very
293:                 ! high and the gradients to very small, so the step is315:                 ! high and the gradients to very small, so the step is
294:                 ! 'converged' but gets rejected.316:                 ! 'converged' but gets rejected.
295:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN317:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN
296:                     CALL REJECT(ENERGY, GRAD)318:                     CALL REJECT(ENERGY, GRAD)
297:                 END IF319:                 END IF
298: 320: 
299:                 ! Reject steps where the cell range has got bigger than the 321:                 ! Reject steps where the cell range has got bigger than the 
871:                         END IF ! GTEST893:                         END IF ! GTEST
872:                     END IF ! End if less than cutoff894:                     END IF ! End if less than cutoff
873:                 END DO ! z loop895:                 END DO ! z loop
874:             END DO ! y loop896:             END DO ! y loop
875:         END DO ! x loop897:         END DO ! x loop
876: 898: 
877:     END SUBROUTINE LJ_GAUSS_PER899:     END SUBROUTINE LJ_GAUSS_PER
878: 900: 
879: !-------------------------------------------------------------------------------901: !-------------------------------------------------------------------------------
880: !902: !
 903: ! Calculates factors for the lattice derivatives which only depend on the cell
 904: ! parameters.
 905: ! ANGLES: alpha, beta and gamma lattice parameters
 906: ! SIZES: a, b and c lattice parameters
 907: !
 908: !-------------------------------------------------------------------------------
 909: 
 910:     PURE TYPE(BOX_DERIV) FUNCTION CALC_BOX_DERIV(ANGLES, SIZES)
 911: 
 912:         IMPLICIT NONE
 913: 
 914:         DOUBLE PRECISION, INTENT(IN) :: ANGLES(3), SIZES(3)
 915: 
 916:         ! V: factor, related to (but not quite equal to) the unit cell volume
 917:         ! C, S: cos and sin of the lattice angles
 918:         ! ORTHOG, DERIV: parts of the output
 919:         DOUBLE PRECISION :: V, C(3), S(3)
 920:         DOUBLE PRECISION :: ORTHOG(3, 3), DERIV(3, 3, 6)
 921: 
 922:         ! Initialise output variables
 923:         ORTHOG(:,:) = 0.D0
 924:         DERIV(:,:,:) = 0.D0
 925: 
 926:         ! Calculate factors
 927:         C(:) = COS(ANGLES(:))
 928:         S(:) = SIN(ANGLES(:))
 929:         V = SQRT(1 - C(1)**2 - C(2)**2 - C(3)**2 + 2.D0 * C(1) * C(2) * C(3))
 930: 
 931:         ! Calculate the orthogonal transformation matrix
 932:         ORTHOG(1, 1) = SIZES(1)
 933:         ORTHOG(1, 2) = SIZES(2) * C(3)
 934:         ORTHOG(1, 3) = SIZES(3) * C(2)
 935:         ORTHOG(2, 2) = SIZES(2) * S(3)
 936:         ORTHOG(2, 3) = SIZES(3) * (C(1) - C(2) * C(3)) / S(3)
 937:         ORTHOG(3, 3) = SIZES(3) * V / S(3)
 938: 
 939:         ! Calculate the derivatives of the othogonal matrix
 940:         ! Derivatives of the top row wrt angles
 941:         DERIV(1, 3, 2) = - SIZES(3) * S(2)
 942:         DERIV(1, 2, 3) = - SIZES(2) * S(3)
 943:         ! Derivatives of the top row wrt to lengths
 944:         DERIV(1, 1, 4) = 1.D0
 945:         DERIV(1, 2, 5) = C(3)
 946:         DERIV(1, 3, 6) = C(2)
 947:         ! Derivatives of the middle row wrt angles
 948:         DERIV(2, 3, 1) = - SIZES(3) * S(1) / S(3)
 949:         DERIV(2, 3, 2) = SIZES(3) * S(2) * C(3) / S(3)
 950:         DERIV(2, 2, 3) = SIZES(2) * C(3)
 951:         DERIV(2, 3, 3) = SIZES(3) * (C(2) - C(1) * C(3)) / S(3)**2
 952:         ! Derivatives of the middle row wrt to lengths
 953:         DERIV(2, 2, 5) = S(3)
 954:         DERIV(2, 3, 6) = (C(1) - C(2) * C(3)) / S(3)
 955:         ! Derivatives of the bottom row wrt angles
 956:         DERIV(3, 3, 1) = SIZES(3) * S(1) * (C(1) - C(2) * C(3)) / (S(3) * V)
 957:         DERIV(3, 3, 2) = SIZES(3) * S(2) * (C(2) - C(3) * C(1)) / (S(3) * V)
 958:         DERIV(3, 3, 3) = SIZES(3) * ((C(3) - C(1) * C(2))/V - C(3) * V/S(3)**2)
 959:         ! Derivatives of the bottom row wrt lengths
 960:         DERIV(3, 3, 6) = V / S(3)
 961: 
 962:         ! Calculate the derivative of wrt the angles
 963:         CALC_BOX_DERIV%V_DERIV(1) = S(1) * (C(1) - C(2) * C(3)) / V
 964:         CALC_BOX_DERIV%V_DERIV(2) = S(2) * (C(2) - C(3) * C(1)) / V
 965:         CALC_BOX_DERIV%V_DERIV(3) = S(3) * (C(3) - C(1) * C(2)) / V
 966: 
 967:         ! Calculate the lengths of the reciprocal lattice vectors
 968:         CALC_BOX_DERIV%RECIP(:) = S(:) / (SIZES(:) * V)
 969: 
 970:         ! Calculate the number of cells to check in each direction
 971:         CALC_BOX_DERIV%CELL_RANGE(:) = CEILING(LJ_GAUSS_RCUT * &
 972:             CALC_BOX_DERIV%RECIP(:))
 973: 
 974:         ! Set the output
 975:         CALC_BOX_DERIV%ORTHOG(:,:) = ORTHOG(:,:)
 976:         CALC_BOX_DERIV%DERIV(:,:,:) = DERIV(:,:,:)
 977:         CALC_BOX_DERIV%V = V
 978: 
 979:     END FUNCTION CALC_BOX_DERIV
 980: 
 981: !-------------------------------------------------------------------------------
 982: !
881: ! Checks whether the unit cell angles are allowed. It is possible to generate983: ! Checks whether the unit cell angles are allowed. It is possible to generate
882: ! sets of angles which do not correspond to a physically possible unit cell.984: ! sets of angles which do not correspond to a physically possible unit cell.
883: ! This would manifest by an imaginary unit cell volume. See985: ! This would manifest by an imaginary unit cell volume. See
884: ! https://journals.iucr.org/a/issues/2011/01/00/au5114/au5114.pdf986: ! https://journals.iucr.org/a/issues/2011/01/00/au5114/au5114.pdf
885: ! 'On the allowed values for the triclinic unit-cell angles',987: ! 'On the allowed values for the triclinic unit-cell angles',
886: ! J. Foadi and G. Evans, Acta Cryst. (2011) A67, 93–95988: ! J. Foadi and G. Evans, Acta Cryst. (2011) A67, 93–95
887: ! ANGLES: the three unit cell angles989: ! ANGLES: the three unit cell angles
888: !990: !
889: !-------------------------------------------------------------------------------991: !-------------------------------------------------------------------------------
890: 992: 
926: ! GRAD_PARAMS: optional, derivatives of potential parameters1028: ! GRAD_PARAMS: optional, derivatives of potential parameters
927: !1029: !
928: !-------------------------------------------------------------------------------1030: !-------------------------------------------------------------------------------
929: 1031: 
930:     SUBROUTINE LJ_GAUSS_PER_TRI(IDI, IDJ, COORDS, GRAD_I, GRAD_J, ENERGY, &1032:     SUBROUTINE LJ_GAUSS_PER_TRI(IDI, IDJ, COORDS, GRAD_I, GRAD_J, ENERGY, &
931:                                 GTEST, GRAD_BOX, BD, GRAD_PARAMS)1033:                                 GTEST, GRAD_BOX, BD, GRAD_PARAMS)
932: 1034: 
933:         ! NATOMS: number of coordinates (including lattice parameters)1035:         ! NATOMS: number of coordinates (including lattice parameters)
934:         USE COMMONS, ONLY: NATOMS1036:         USE COMMONS, ONLY: NATOMS
935: 1037: 
936:         ! BOX_DERIV: type for storing box derivative information 
937:         USE TRICLINIC_MOD, ONLY: BOX_DERIV 
938:  
939:         IMPLICIT NONE1038:         IMPLICIT NONE
940: 1039: 
941:         INTEGER, INTENT(IN) :: IDI, IDJ1040:         INTEGER, INTENT(IN) :: IDI, IDJ
942:         DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS)1041:         DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS)
943:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_I(3), GRAD_J(3), ENERGY1042:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_I(3), GRAD_J(3), ENERGY
944:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_BOX(6)1043:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_BOX(6)
945:         LOGICAL, INTENT(IN) :: GTEST1044:         LOGICAL, INTENT(IN) :: GTEST
946:         TYPE(BOX_DERIV), INTENT(IN) :: BD1045:         TYPE(BOX_DERIV), INTENT(IN) :: BD
947:         DOUBLE PRECISION, INTENT(INOUT), OPTIONAL :: GRAD_PARAMS(3)1046:         DOUBLE PRECISION, INTENT(INOUT), OPTIONAL :: GRAD_PARAMS(3)
948: 1047: 
1130: ! adopting physically impossible angle combinations. Uses a WCA potential1229: ! adopting physically impossible angle combinations. Uses a WCA potential
1131: ! ENERGY: energy of the system1230: ! ENERGY: energy of the system
1132: ! GRAD_ANGLES: gradients for the box angles1231: ! GRAD_ANGLES: gradients for the box angles
1133: ! BD: information on the box and its derivatives1232: ! BD: information on the box and its derivatives
1134: ! GTEST: whether to calculate gradients1233: ! GTEST: whether to calculate gradients
1135: !1234: !
1136: !-------------------------------------------------------------------------------1235: !-------------------------------------------------------------------------------
1137: 1236: 
1138:     SUBROUTINE CONSTRAIN_VOLUME(ENERGY, GRAD_ANGLES, BD, GTEST)1237:     SUBROUTINE CONSTRAIN_VOLUME(ENERGY, GRAD_ANGLES, BD, GTEST)
1139: 1238: 
1140:         ! BOX_DERIV: type for storing box derivative information 
1141:         USE TRICLINIC_MOD, ONLY: BOX_DERIV 
1142:  
1143:         IMPLICIT NONE1239:         IMPLICIT NONE
1144: 1240: 
1145:         DOUBLE PRECISION, INTENT(INOUT) :: ENERGY, GRAD_ANGLES(3)1241:         DOUBLE PRECISION, INTENT(INOUT) :: ENERGY, GRAD_ANGLES(3)
1146:         TYPE(BOX_DERIV), INTENT(IN) :: BD1242:         TYPE(BOX_DERIV), INTENT(IN) :: BD
1147:         LOGICAL, INTENT(IN) :: GTEST1243:         LOGICAL, INTENT(IN) :: GTEST
1148: 1244: 
 1245: !        WRITE (*, *) 'BD%V = ', BD%V
 1246: !        WRITE (*, *) 'BD%V_DERIV(:) = ', BD%V_DERIV(:)
 1247: 
1149:         ! Add a purely repulsive (away from zero) WCA energy term1248:         ! Add a purely repulsive (away from zero) WCA energy term
1150:         IF (BD%V .LT. V_SIGMA * 2.D0**(1.D0/6.D0)) THEN1249:         IF (BD%V .LT. V_SIGMA**(1.D0/6.D0)) THEN
1151:             ENERGY = ENERGY + 4.D0 * V_EPS * &1250:             ENERGY = ENERGY + 4.D0 * V_EPS * &
1152:                      ((V_SIGMA / BD%V)**(12) - (V_SIGMA / BD%V)**(6)) + V_EPS1251:                      ((V_SIGMA / BD%V)**(12) - (V_SIGMA / BD%V)**(6)) + V_EPS
1153: 1252: 
1154:             IF (GTEST) THEN1253:             IF (GTEST) THEN
1155:                 ! Add the gradients, if necessary1254:                 ! Add the gradients, if necessary
1156:                 GRAD_ANGLES(:) = GRAD_ANGLES(:) + 24.D0 * V_EPS / V_SIGMA * &1255:                 GRAD_ANGLES(:) = GRAD_ANGLES(:) + 24.D0 * V_EPS / V_SIGMA * &
1157:                     ((V_SIGMA / BD%V)**(7) - 2.D0 * (V_SIGMA / BD%V)**(13)) * &1256:                     ((V_SIGMA / BD%V)**(7) - 2.D0 * (V_SIGMA / BD%V)**(13)) * &
1158:                     BD%V_DERIV(:)1257:                     BD%V_DERIV(:)
1159:             END IF1258:             END IF
1160:         END IF1259:         END IF
1252:         ! NSAVE: number of minima to write1351:         ! NSAVE: number of minima to write
1253:         ! BOXx: Box dimensions1352:         ! BOXx: Box dimensions
1254:         ! PERIODIC: whether periodic boundary conditions are in use1353:         ! PERIODIC: whether periodic boundary conditions are in use
1255:         USE COMMONS, ONLY: NATOMS, NSAVE, PERIODIC, BOXLX, BOXLY, BOXLZ1354:         USE COMMONS, ONLY: NATOMS, NSAVE, PERIODIC, BOXLX, BOXLY, BOXLZ
1256: 1355: 
1257:         ! QMIN: energies of the saved minima1356:         ! QMIN: energies of the saved minima
1258:         ! QMINP: coordinates of the saved minima1357:         ! QMINP: coordinates of the saved minima
1259:         ! FF: step at which the saved minima was first found1358:         ! FF: step at which the saved minima was first found
1260:         USE QMODULE, ONLY: QMIN, QMINP, FF1359:         USE QMODULE, ONLY: QMIN, QMINP, FF
1261: 1360: 
1262:         ! BOX_DERIV: type for storing box derivative information 
1263:         ! CALC_BOX_DERIV: calculates box derivative information 
1264:         USE TRICLINIC_MOD, ONLY: BOX_DERIV, CALC_BOX_DERIV 
1265:  
1266:         IMPLICIT NONE1361:         IMPLICIT NONE
1267: 1362: 
1268:         ! GETUNIT: function for fiding a free file unit1363:         ! GETUNIT: function for fiding a free file unit
1269:         ! FUNIT: output file unit1364:         ! FUNIT: output file unit
1270:         ! NMOL: actual number of particles1365:         ! NMOL: actual number of particles
1271:         ! I, J: loop iterators1366:         ! I, J: loop iterators
1272:         INTEGER :: GETUNIT, FUNIT, NMOL, I, J1367:         INTEGER :: GETUNIT, FUNIT, NMOL, I, J
1273: 1368: 
1274:         ! BOX_SIZE: convenience array for the box sizes1369:         ! BOX_SIZE: convenience array for the box sizes
1275:         ! BOX_ANGLES: convenience array for the box angles1370:         ! BOX_ANGLES: convenience array for the box angles


r33077/opp.f90 2017-07-26 19:30:28.920733101 +0100 r33076/opp.f90 2017-07-26 19:30:31.120762315 +0100
 24: ! separated out and generalised to all isotropic pair potentials. 24: ! separated out and generalised to all isotropic pair potentials.
 25: ! Questions to jwrm2. 25: ! Questions to jwrm2.
 26: ! 26: !
 27: MODULE OPP_MOD 27: MODULE OPP_MOD
 28:  28: 
 29: ! Public parameters 29: ! Public parameters
 30: PUBLIC :: OPP_MODE, OPP_RCUT, OPP_K, OPP_PHI, OPP_PARAMS 30: PUBLIC :: OPP_MODE, OPP_RCUT, OPP_K, OPP_PHI, OPP_PARAMS
 31: ! Public functions 31: ! Public functions
 32: PUBLIC :: OPP, OPP_INITIALISE, VIEW_OPP, OPP_TAKESTEP 32: PUBLIC :: OPP, OPP_INITIALISE, VIEW_OPP, OPP_TAKESTEP
 33: ! Private parameters 33: ! Private parameters
 34: PRIVATE :: V_EPS, V_SIGMA, FD_RCUT, IMAGE_CUTOFF!, SD_RCUT 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 35: PRIVATE :: MAX_LENGTH_STEP, MAX_ANGLE_STEP, MAX_K_STEP, MAX_PHI_STEP
 36: PRIVATE :: K_REPULSION, K_LOWER, K_UPPER 36: PRIVATE :: K_REPULSION, K_LOWER, K_UPPER
 37: ! Private functions 37: ! Private functions
 38: PRIVATE :: PERIODIC_RESET, OPP_PER, OPP_PER_TRI, CALC_ENERGY 38: PRIVATE :: PERIODIC_RESET, OPP_PER, OPP_PER_TRI, CALC_ENERGY
 39: PRIVATE :: CALC_GRAD, CALC_SEC, CALC_FCTS, CONSTRAIN_PARAMETERS 39: PRIVATE :: CALC_GRAD, CALC_FCTS, CONSTRAIN_PARAMETERS
 40: PRIVATE :: CHECK_ANGLES, REJECT, CONSTRAIN_VOLUME 40: PRIVATE :: BOX_DERIV, CALC_BOX_DERIV, CHECK_ANGLES, REJECT, CONSTRAIN_VOLUME
 41:  41: 
 42: ! K: frequency parameter, controls the number of wells 42: ! K: frequency parameter, controls the number of wells
 43: ! PHI: phase parameter, controls the position of the first minimum 43: ! PHI: phase parameter, controls the position of the first minimum
 44: ! RCUT: cutoff distance at which the energy, gradients and Hessian of 44: ! RCUT: cutoff distance at which the energy, gradients and Hessian of
 45: !       the potential go smoothly to zero. 45: !       the potential go smoothly to zero.
 46: DOUBLE PRECISION :: OPP_K, OPP_PHI, OPP_RCUT 46: DOUBLE PRECISION :: OPP_K, OPP_PHI, OPP_RCUT
 47:  47: 
 48: ! OPP_MODE: 0 standard minimisation 48: ! OPP_MODE: 0 standard minimisation
 49: !           1 the final triplet of coordinates are the unit cell length 49: !           1 the final triplet of coordinates are the unit cell length
 50: !             parameters and the second last triplet are the unit cell 50: !             parameters and the second last triplet are the unit cell
 89: ! V_EPS: scaling factor for the unit cell volume contraint energy 89: ! V_EPS: scaling factor for the unit cell volume contraint energy
 90: ! V_SIGMA: WCA cutoff distance 90: ! V_SIGMA: WCA cutoff distance
 91: ! K_LOWER, K_UPPER: bounds on OPP_K 91: ! K_LOWER, K_UPPER: bounds on OPP_K
 92: ! K_REPULSION: harmonic force constant for 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 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 94: DOUBLE PRECISION, PARAMETER :: V_EPS = 1.D-3
 95: DOUBLE PRECISION, PARAMETER :: V_SIGMA = 3.D-1 95: DOUBLE PRECISION, PARAMETER :: V_SIGMA = 3.D-1
 96: DOUBLE PRECISION, PARAMETER :: K_LOWER = 5.D0, K_UPPER = 20.D0 96: DOUBLE PRECISION, PARAMETER :: K_LOWER = 5.D0, K_UPPER = 20.D0
 97: DOUBLE PRECISION, PARAMETER :: K_REPULSION = 1.D4 97: DOUBLE PRECISION, PARAMETER :: K_REPULSION = 1.D4
 98:  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: 
 99: CONTAINS125: CONTAINS
100: 126: 
101: !-------------------------------------------------------------------------------127: !-------------------------------------------------------------------------------
102: !128: !
103: ! Main routine for calculating the energy and gradients.129: ! Main routine for calculating the energy and gradients.
104: ! X: coordinates array130: ! X: coordinates array
105: ! GRAD: gradients array131: ! GRAD: gradients array
106: ! ENERGY: calculated energy132: ! ENERGY: calculated energy
107: ! GTEST: whether gradients should be calculated133: ! GTEST: whether gradients should be calculated
108: ! STEST: whether the Hessian should be calculated 
109: !134: !
110: !-------------------------------------------------------------------------------135: !-------------------------------------------------------------------------------
111:     SUBROUTINE OPP(X, GRAD, ENERGY, GTEST, STEST)136:     SUBROUTINE OPP(X, GRAD, ENERGY, GTEST)
112: 137: 
113:         ! MYUNIT: file unit for main output file138:         ! MYUNIT: file unit for main output file
114:         ! NATOMS: number of particles139:         ! NATOMS: number of particles
115:         ! PERIODIC: whether to use periodic boundary conditions140:         ! PERIODIC: whether to use periodic boundary conditions
116:         USE COMMONS, ONLY: MYUNIT, NATOMS, PERIODIC141:         USE COMMONS, ONLY: MYUNIT, NATOMS, PERIODIC
117: 142: 
118:         ! HESS; the Hessian matrix 
119:         USE MODHESS, ONLY: HESS 
120:  
121:         ! BOX_DERIV: type for storing box derivative information 
122:         ! CALC_BOX_DERIV: calculates box derivative information 
123:         USE TRICLINIC_MOD, ONLY: BOX_DERIV, CALC_BOX_DERIV 
124:  
125:         IMPLICIT NONE143:         IMPLICIT NONE
126: 144: 
127:         DOUBLE PRECISION, INTENT(INOUT) :: X(3*NATOMS)145:         DOUBLE PRECISION, INTENT(INOUT) :: X(3*NATOMS)
128:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOMS), ENERGY146:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOMS), ENERGY
129:         LOGICAL, INTENT(IN) :: GTEST, STEST147:         LOGICAL, INTENT(IN) :: GTEST
130: 148: 
131:         ! I, J: indices for looping over particles149:         ! I, J: indices for looping over particles
132:         ! NMOL: actual number of particles150:         ! NMOL: actual number of particles
133:         INTEGER :: I, J, NMOL151:         INTEGER :: I, J, NMOL
134: 152: 
135:         ! RIJ: Interparticle vector153:         ! RIJ: Interparticle vector
136:         ! MODRIJ: distance between a pair of particles154:         ! MODRIJ: distance between a pair of particles
137:         ! DVDR: derivative of pair potential wrt MODRIJ155:         ! DVDR: derivative of pair potential wrt MODRIJ
138:         ! D2VDR2: D(DVDR/MODRIJ)/DMODRIJ * 1/MODRIJ156:         DOUBLE PRECISION :: RIJ(3), MODRIJ, DVDR
139:         DOUBLE PRECISION :: RIJ(3), MODRIJ, DVDR, D2VDR2 
140: 157: 
141:         ! Factors for triclinic lattice158:         ! Factors for triclinic lattice
142:         TYPE(BOX_DERIV) :: BD159:         TYPE(BOX_DERIV) :: BD
143: 160: 
144:         ! Initialise output variables161:         ! Initialise output variables
145:         ENERGY = 0.D0162:         ENERGY = 0.D0
146:         GRAD(:) = 0.D0163:         GRAD(:) = 0.D0
147:         IF (STEST) HESS(:,:) = 0.D0 
148: 164: 
149:         ! Calculate factors that do depend on parameters, but not on RIJ165:         ! Calculate factors that do depend on parameters, but not on RIJ
150:         IF (OPP_MODE .EQ. 2) THEN166:         IF (OPP_MODE .EQ. 2) THEN
151:             CALL CALC_FCTS(X(NATOMS*3-8:NATOMS*3-6))167:             CALL CALC_FCTS(X(NATOMS*3-8:NATOMS*3-6))
152:         ELSE168:         ELSE
153:             CALL CALC_FCTS()169:             CALL CALC_FCTS()
154:         END IF170:         END IF
155: 171: 
156:         IF (.NOT. PERIODIC) THEN172:         IF (.NOT. PERIODIC) THEN
157:             IF (OPP_MODE .EQ. 0) THEN173:             IF (OPP_MODE .EQ. 0) THEN
158:                 ! Normal cluster174:                 ! Normal cluster
159:                 NMOL = NATOMS175:                 NMOL = NATOMS
160:                 DO I = 1, NMOL ! Outer loop over particles176:                 DO I = 1, NMOL-1 ! Outer loop over particles
161:                     DO J = 1, NMOL ! Inner loop over particles177:                     DO J = I+1, NMOL ! Inner loop over particles
162:  
163:                         ! No self interaction 
164:                         IF (I .EQ. J) CYCLE 
165:  
166:                         ! Get the particle-particle distance178:                         ! Get the particle-particle distance
167:                         RIJ(:) = X(3*I-2:3*I) - X(3*J-2:3*J)179:                         RIJ(:) = X(3*I-2:3*I) - X(3*J-2:3*J)
168:                         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))180:                         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))
169: 181: 
170:                         ! Add energy and gradients182:                         ! Add energy and gradients
171:                         IF (MODRIJ .LT. OPP_RCUT) THEN183:                         IF (MODRIJ .LT. OPP_RCUT) THEN
172: 184:                             ENERGY = ENERGY + CALC_ENERGY(MODRIJ)
173:                             ! Don't double count energy 
174:                             IF (J .GT. I) THEN 
175:                                 ENERGY = ENERGY + CALC_ENERGY(MODRIJ) 
176:                             END IF 
177: 185: 
178:                             IF (GTEST) THEN186:                             IF (GTEST) THEN
179:                                 DVDR = CALC_GRAD(MODRIJ) / MODRIJ187:                                 DVDR = CALC_GRAD(MODRIJ)
180:                                 GRAD(3*I-2:3*I) = GRAD(3*I-2:3*I) + DVDR*RIJ(:)188:                                 GRAD(3*I-2:3*I) = GRAD(3*I-2:3*I) + &
181: !                                GRAD(3*J-2:3*J) = GRAD(3*J-2:3*J) - DVDR*RIJ(:)189:                                                   DVDR*RIJ(:)/MODRIJ
 190:                                 GRAD(3*J-2:3*J) = GRAD(3*J-2:3*J) - &
 191:                                                   DVDR*RIJ(:)/MODRIJ
182:                             END IF ! GTEST192:                             END IF ! GTEST
183:  
184:                             IF (STEST) THEN 
185:                                 ! Hessian terms 
186:                                 DVDR = DVDR / MODRIJ 
187:                                 D2VDR2 = (CALC_SEC(MODRIJ) - DVDR) / MODRIJ**2 
188:                                 ! Same atom, same coordinate 
189:                                 HESS(3*I-2,3*I-2) = HESS(3*I-2,3*I-2) + & 
190:                                     D2VDR2 * RIJ(1)**2 + DVDR 
191:                                 HESS(3*I-1,3*I-1) = HESS(3*I-1,3*I-1) + & 
192:                                     D2VDR2 * RIJ(2)**2 + DVDR 
193:                                 HESS(3*I  ,3*I  ) = HESS(3*I  ,3*I  ) + & 
194:                                     D2VDR2 * RIJ(3)**2 + DVDR 
195:                                 ! Same atom, different coordinate 
196:                                 HESS(3*I-2,3*I-1) = HESS(3*I-2,3*I-1) + & 
197:                                     D2VDR2 * RIJ(1) * RIJ(2) 
198:                                 HESS(3*I-1,3*I  ) = HESS(3*I-1,3*I  ) + & 
199:                                     D2VDR2 * RIJ(2) * RIJ(3) 
200:                                 HESS(3*I  ,3*I-2) = HESS(3*I  ,3*I-2) + & 
201:                                     D2VDR2 * RIJ(3) * RIJ(1) 
202:                                 HESS(3*I-1,3*I-2) = HESS(3*I-1,3*I-2) + & 
203:                                     D2VDR2 * RIJ(2) * RIJ(1) 
204:                                 HESS(3*I  ,3*I-1) = HESS(3*I  ,3*I-1) + & 
205:                                     D2VDR2 * RIJ(3) * RIJ(2) 
206:                                 HESS(3*I-2,3*I  ) = HESS(3*I-2,3*I) + & 
207:                                     D2VDR2 * RIJ(1) * RIJ(3) 
208:                                 ! Different atom, same coordinate 
209:                                 HESS(3*I-2,3*J-2) = - D2VDR2 * RIJ(1)**2 - DVDR 
210:                                 HESS(3*I-1,3*J-1) = - D2VDR2 * RIJ(2)**2 - DVDR 
211:                                 HESS(3*I  ,3*J  ) = - D2VDR2 * RIJ(3)**2 - DVDR 
212:                                 HESS(3*J-2,3*I-2) = - D2VDR2 * RIJ(1)**2 - DVDR 
213:                                 HESS(3*J-1,3*I-1) = - D2VDR2 * RIJ(2)**2 - DVDR 
214:                                 HESS(3*J  ,3*I  ) = - D2VDR2 * RIJ(3)**2 - DVDR 
215:                                 ! Different atom, different coordinate 
216:                                 HESS(3*I-2,3*J-1) = - D2VDR2 * RIJ(1) * RIJ(2) 
217:                                 HESS(3*I-1,3*J  ) = - D2VDR2 * RIJ(2) * RIJ(3) 
218:                                 HESS(3*I-2,3*J  ) = - D2VDR2 * RIJ(3) * RIJ(1) 
219:                                 HESS(3*J-1,3*I-2) = - D2VDR2 * RIJ(1) * RIJ(2) 
220:                                 HESS(3*J  ,3*I-1) = - D2VDR2 * RIJ(2) * RIJ(3) 
221:                                 HESS(3*J  ,3*I-2) = - D2VDR2 * RIJ(3) * RIJ(1) 
222:                             END IF ! End second derivative 
223:                         END IF ! If less than cutoff distance193:                         END IF ! If less than cutoff distance
224:                     END DO ! Inner loop over particles194:                     END DO ! Inner loop over particles
225:                 END DO ! Outer loop over particles195:                 END DO ! Outer loop over particles
226:             ELSE ! mode not equal to zero196:             ELSE ! mode not equal to zero
227:                 WRITE (MYUNIT, *) 'OPP> ERROR, PERIODIC must be ', &197:                 WRITE (MYUNIT, *) 'OPP> ERROR, PERIODIC must be ', &
228:                                   'specified if mode is not 0'198:                                   'specified if mode is not 0'
229:                 STOP199:                 STOP
230:             END IF200:             END IF
231:         ELSE ! PERIODIC201:         ELSE ! PERIODIC
232:             SELECT CASE (OPP_MODE)202:             SELECT CASE (OPP_MODE)
233:             CASE(0)203:             CASE(0)
234:                 ! Periodic, fixed box204:                 ! Periodic, fixed box
235:                 ! Reset all particles to within the box205:                 ! Reset all particles to within the box
236:                 CALL PERIODIC_RESET(X(:))206:                 CALL PERIODIC_RESET(X(:))
237: 207: 
238:                 NMOL = NATOMS208:                 NMOL = NATOMS
239:                 ! We need to include self-interactions of particles in different209:                 ! We need to include self-interactions of particles in different
240:                 ! unit cells210:                 ! unit cells
241:                 DO I = 1, NMOL ! Outer loop over particles211:                 DO I = 1, NMOL ! Outer loop over particles
242:                     DO J = 1, NMOL ! Inner loop over particles212:                     DO J = I, NMOL ! Inner loop over particles
243:                         ! Add energy and gradients213:                         ! Add energy and gradients
244:                         CALL OPP_PER(I, J, X(3*I-2:3*I), X(3*J-2:3*J), &214:                         CALL OPP_PER(I, J, X(3*I-2:3*I), X(3*J-2:3*J), &
245:                                      GRAD(3*I-2:3*I), ENERGY, GTEST, STEST)215:                                      GRAD(3*I-2:3*I), GRAD(3*J-2:3*J), &
 216:                                      ENERGY, GTEST)
246:                     END DO ! Inner loop over particles217:                     END DO ! Inner loop over particles
247:                 END DO ! Outer loop over particles218:                 END DO ! Outer loop over particles
248: 219: 
249:             CASE(1)220:             CASE(1)
250:                 ! Triclinic varying lattice221:                 ! Triclinic varying lattice
251:                 ! Reset all particles to within the box222:                 ! Reset all particles to within the box
252:                 CALL PERIODIC_RESET(X(:))223:                 CALL PERIODIC_RESET(X(:))
253: 224: 
254: !                DO I = 1, NATOMS225: !                DO I = 1, NATOMS
255: !                    WRITE (MYUNIT, *) X(3*I-2:3*I)226: !                    WRITE (MYUNIT, *) X(3*I-2:3*I)
256: !               END DO227: !                END DO
257: 228: 
258:                 ! Calculate box derivative factors229:                 ! Calculate box derivative factors
259:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &230:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &
260:                                     X(3*NATOMS-2:3*NATOMS  ), OPP_RCUT)231:                                     X(3*NATOMS-2:3*NATOMS  ))
261: 232: 
262:                 ! Check whether the unit cell angles are physically possible.233:                 ! Check whether the unit cell angles are physically possible.
263:                 ! Reject if not.234:                 ! Reject if not.
264:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN235:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN
265:                     CALL REJECT(ENERGY, GRAD)236:                     CALL REJECT(ENERGY, GRAD)
266:                     RETURN237:                     RETURN
267:                 END IF238:                 END IF
268: 239: 
269:                 ! Reject steps where the cell range has got bigger than the 240:                 ! Reject steps where the cell range has got bigger than the 
270:                 ! cutoff. Such unit cells are highly distorted and there are 241:                 ! cutoff. Such unit cells are highly distorted and there are 
275:                     RETURN246:                     RETURN
276:                 END IF247:                 END IF
277: 248: 
278:                 NMOL = NATOMS-2249:                 NMOL = NATOMS-2
279:                 ! We need to include self-interactions of particles in different250:                 ! We need to include self-interactions of particles in different
280:                 ! unit cells251:                 ! unit cells
281:                 DO I = 1, NMOL! Outer loop over particles252:                 DO I = 1, NMOL! Outer loop over particles
282:                     DO J = I, NMOL ! Inner loop over particles253:                     DO J = I, NMOL ! Inner loop over particles
283:                         ! Add energy and gradients254:                         ! Add energy and gradients
284:                         CALL OPP_PER_TRI(I, J, X(:), GRAD(3*I-2:3*I), &255:                         CALL OPP_PER_TRI(I, J, X(:), GRAD(3*I-2:3*I), &
285:                             GRAD(3*J-2:3*J), ENERGY, GTEST, STEST, &256:                             GRAD(3*J-2:3*J), ENERGY, GTEST, &
286:                             GRAD(3*NATOMS-5:3*NATOMS), BD)257:                             GRAD(3*NATOMS-5:3*NATOMS), BD)
287:                     END DO ! Inner loop over particles258:                     END DO ! Inner loop over particles
288:                 END DO ! Outer loop over particles259:                 END DO ! Outer loop over particles
289: 260: 
290:                 ! Constrain the box volume to be greater than 0 with a WCA-style261:                 ! Constrain the box volume to be greater than 0 with a WCA-style
291:                 ! repulsion262:                 ! repulsion
292:                 CALL CONSTRAIN_VOLUME(ENERGY, GRAD(3*NATOMS-5:3*NATOMS-3), BD, &263:                 CALL CONSTRAIN_VOLUME(ENERGY, GRAD(3*NATOMS-5:3*NATOMS-3), BD, &
293:                                       GTEST, STEST)264:                                       GTEST)
294: 265: 
295:             CASE(2)266:             CASE(2)
296:                 ! Triclinic varying lattice, with varying parameters267:                 ! Triclinic varying lattice, with varying parameters
297:                 ! Reset all particles to within the box268:                 ! Reset all particles to within the box
298:                 CALL PERIODIC_RESET(X(:))269:                 CALL PERIODIC_RESET(X(:))
299: 270: 
300:                 ! Calculate box derivative factors271:                 ! Calculate box derivative factors
301:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &272:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &
302:                                     X(3*NATOMS-2:3*NATOMS  ), OPP_RCUT)273:                                     X(3*NATOMS-2:3*NATOMS  ))
303: 274: 
304:                 ! Check whether the unit cell angles are physically possible.275:                 ! Check whether the unit cell angles are physically possible.
305:                 ! If we have a disallowed unit cell, we set the energy to very276:                 ! If we have a disallowed unit cell, we set the energy to very
306:                 ! high and the gradients to very small, so the step is277:                 ! high and the gradients to very small, so the step is
307:                 ! 'converged' but gets rejected.278:                 ! 'converged' but gets rejected.
308:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN279:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN
309:                     CALL REJECT(ENERGY, GRAD)280:                     CALL REJECT(ENERGY, GRAD)
310:                 END IF281:                 END IF
311: 282: 
312:                 ! Reject steps where the cell range has got bigger than the 283:                 ! Reject steps where the cell range has got bigger than the 
318:                     RETURN289:                     RETURN
319:                 END IF290:                 END IF
320: 291: 
321:                 NMOL = NATOMS-3292:                 NMOL = NATOMS-3
322:                 ! We need to include self-interactions of particles in different293:                 ! We need to include self-interactions of particles in different
323:                 ! unit cells294:                 ! unit cells
324:                 DO I = 1, NMOL! Outer loop over particles295:                 DO I = 1, NMOL! Outer loop over particles
325:                     DO J = I, NMOL ! Inner loop over particles296:                     DO J = I, NMOL ! Inner loop over particles
326:                         ! Add energy and gradients297:                         ! Add energy and gradients
327:                         CALL OPP_PER_TRI(I, J, X(:), GRAD(3*I-2:3*I), &298:                         CALL OPP_PER_TRI(I, J, X(:), GRAD(3*I-2:3*I), &
328:                             GRAD(3*J-2:3*J), ENERGY, GTEST, STEST, &299:                             GRAD(3*J-2:3*J), ENERGY, GTEST, &
329:                             GRAD(3*NATOMS-5:3*NATOMS), BD, &300:                             GRAD(3*NATOMS-5:3*NATOMS), BD, &
330:                             GRAD(3*NATOMS-8:3*NATOMS-6))301:                             GRAD(3*NATOMS-8:3*NATOMS-6))
331:                     END DO ! Inner loop over particles302:                     END DO ! Inner loop over particles
332:                 END DO ! Outer loop over particles303:                 END DO ! Outer loop over particles
333: 304: 
334:                 ! Constrain the box volume to be greater than 0 with a WCA-style305:                 ! Constrain the box volume to be greater than 0 with a WCA-style
335:                 ! repulsion306:                 ! repulsion
336:                 CALL CONSTRAIN_VOLUME(ENERGY, GRAD(3*NATOMS-5:3*NATOMS-3), BD, &307:                 CALL CONSTRAIN_VOLUME(ENERGY, GRAD(3*NATOMS-5:3*NATOMS-3), BD, &
337:                                       GTEST, STEST)308:                                       GTEST)
338: 309: 
339:                 ! Constrain the potential parameters with a harmonic repulsion310:                 ! Constrain the potential parameters with a harmonic repulsion
340:                 CALL CONSTRAIN_PARAMETERS(ENERGY, GRAD(3*NATOMS-8:3*NATOMS-6), &311:                 CALL CONSTRAIN_PARAMETERS(ENERGY, GRAD(3*NATOMS-8:3*NATOMS-6), &
341:                                           GTEST, STEST)312:                                           GTEST)
342:             END SELECT ! Mode selections313:             END SELECT ! Mode selections
343:         END IF ! Periodic314:         END IF ! Periodic
344: 315: 
 316:         FLUSH(MYUNIT)
 317: 
345:     END SUBROUTINE OPP318:     END SUBROUTINE OPP
346: 319: 
347: !-------------------------------------------------------------------------------320: !-------------------------------------------------------------------------------
348: !321: !
349: ! Rejects a step by setting the energy very high and the gradients very low.322: ! Rejects a step by setting the energy very high and the gradients very low.
350: ! This tricks L-BFGS into thinking it's converged, but the high energy of the323: ! This tricks L-BFGS into thinking it's converged, but the high energy of the
351: ! quench will ensure it is rejected.324: ! quench will ensure it is rejected.
352: ! ENERGY: system energy325: ! ENERGY: system energy
353: ! GRAD: system gradients326: ! GRAD: system gradients
354: !327: !
469: 442: 
470:         ! Unmodified gradient443:         ! Unmodified gradient
471:         DVDR = - 15.D0 * MODRIJ**(-16) - &444:         DVDR = - 15.D0 * MODRIJ**(-16) - &
472:                OPP_K * MODRIJ**(-3) * SIN(OPP_K * (MODRIJ - 1) - OPP_PHI) - & 445:                OPP_K * MODRIJ**(-3) * SIN(OPP_K * (MODRIJ - 1) - OPP_PHI) - & 
473:                3.D0 * MODRIJ**(-4) * COS(OPP_K * (MODRIJ - 1) - OPP_PHI)446:                3.D0 * MODRIJ**(-4) * COS(OPP_K * (MODRIJ - 1) - OPP_PHI)
474: 447: 
475:         ! Term to make the first derivative go to zero at OPP_RCUT448:         ! Term to make the first derivative go to zero at OPP_RCUT
476:         DVDR = DVDR - FD_RCUT449:         DVDR = DVDR - FD_RCUT
477: 450: 
478:         ! Term to make the second derivative go to zero at OPP_RCUT451:         ! Term to make the second derivative go to zero at OPP_RCUT
479: !        DVDR = DVDR - (MODRIJ - OPP_RCUT) * SD_RCUT452:         DVDR = DVDR - (MODRIJ - OPP_RCUT) * SD_RCUT
480: 453: 
481:         CALC_GRAD = DVDR454:         CALC_GRAD = DVDR
482: 455: 
483:     END FUNCTION CALC_GRAD456:     END FUNCTION CALC_GRAD
484: 457: 
485: !-------------------------------------------------------------------------------458: !-------------------------------------------------------------------------------
486: !459: !
487: ! Calculates the pair potential second derivative for a pair of particles. 
488: ! MODRIJ: distance beteen the particle pair 
489: ! 
490: !------------------------------------------------------------------------------- 
491:     PURE DOUBLE PRECISION FUNCTION CALC_SEC (MODRIJ) 
492:  
493:         IMPLICIT NONE 
494:  
495:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
496:  
497:         ! Unmodified gradient 
498:         CALC_SEC = 240.D0 * MODRIJ**(-17) - & 
499:             OPP_K**2 * MODRIJ**(-3) * COS(OPP_K * (MODRIJ - 1) - OPP_PHI) + &  
500:             6.D0 * OPP_K * MODRIJ**(-4) * SIN(OPP_K * (MODRIJ - 1) - OPP_PHI) +& 
501:             12.D0 * MODRIJ**(-5) * COS(OPP_K * (MODRIJ - 1) - OPP_PHI) 
502:  
503:     END FUNCTION CALC_SEC 
504:  
505: !------------------------------------------------------------------------------- 
506: ! 
507: ! Initialisation. Check values for OPP_MODE and OPP_PARAMS.460: ! Initialisation. Check values for OPP_MODE and OPP_PARAMS.
508: !461: !
509: !-------------------------------------------------------------------------------462: !-------------------------------------------------------------------------------
510:     SUBROUTINE OPP_INITIALISE()463:     SUBROUTINE OPP_INITIALISE()
511: 464: 
512:         ! MYUNIT: file unit for main GMIN output465:         ! MYUNIT: file unit for main GMIN output
513:         USE COMMONS, ONLY: MYUNIT466:         USE COMMONS, ONLY: MYUNIT
514: 467: 
515:         IMPLICIT NONE468:         IMPLICIT NONE
516: 469: 
705: !-------------------------------------------------------------------------------658: !-------------------------------------------------------------------------------
706: !659: !
707: ! Finds the energy and gradients for the othogonal fixed periodic system660: ! Finds the energy and gradients for the othogonal fixed periodic system
708: ! IDI, IDJ: id numbers of the particles661: ! IDI, IDJ: id numbers of the particles
709: ! COORDS_I: coordinates of particle I662: ! COORDS_I: coordinates of particle I
710: ! COORDS_J: coordinates of particle J663: ! COORDS_J: coordinates of particle J
711: ! GRAD_I: gradients for particle I664: ! GRAD_I: gradients for particle I
712: ! GRAD_J: gradients for particle J665: ! GRAD_J: gradients for particle J
713: ! ENERGY: energy contribution for this pair666: ! ENERGY: energy contribution for this pair
714: ! GTEST: whether to calculate gradients667: ! GTEST: whether to calculate gradients
715: ! STEST: whether to calculate Hessian 
716: !668: !
717: !-------------------------------------------------------------------------------669: !-------------------------------------------------------------------------------
718: 670: 
719:     SUBROUTINE OPP_PER(IDI, IDJ, COORDS_I, COORDS_J, GRAD_I, ENERGY, &671:     SUBROUTINE OPP_PER(IDI, IDJ, COORDS_I,COORDS_J,GRAD_I,GRAD_J,ENERGY, GTEST)
720:                        GTEST, STEST) 
721: 672: 
722:         ! BOXLX, BOXLY, BOXLZ: dimensions of the box673:         ! BOXLX, BOXLY, BOXLZ: dimensions of the box
723:         USE COMMONS, ONLY: BOXLX, BOXLY, BOXLZ674:         USE COMMONS, ONLY: BOXLX, BOXLY, BOXLZ
724: 675: 
725:         ! HESS; the Hessian matrix 
726:         USE MODHESS, ONLY: HESS 
727:  
728:         IMPLICIT NONE676:         IMPLICIT NONE
729: 677: 
730:         INTEGER, INTENT(IN) :: IDI, IDJ678:         INTEGER, INTENT(IN) :: IDI, IDJ
731:         DOUBLE PRECISION, INTENT(IN) :: COORDS_I(3), COORDS_J(3)679:         DOUBLE PRECISION, INTENT(IN) :: COORDS_I(3), COORDS_J(3)
732:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_I(3), ENERGY680:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_I(3), GRAD_J(3), ENERGY
733:         LOGICAL, INTENT(IN) :: GTEST, STEST681:         LOGICAL, INTENT(IN) :: GTEST
734: 682: 
735:         ! XJ: copy of j coordinates683:         ! XJ: copy of j coordinates
736:         ! RIJ: particle-particle vector684:         ! RIJ: particle-particle vector
737:         ! MODRIJ: particle-particle distance685:         ! MODRIJ: particle-particle distance
738:         ! DVDR: derivative of pair potential wrt MODRIJ686:         ! DVDR: derivative of pair potential wrt MODRIJ
739:         ! D2VDR2: D(DVDR/MODRIJ)/DMODRIJ * 1/MODRIJ 
740:         ! BOX_SIZE: convenience copy of the box size687:         ! BOX_SIZE: convenience copy of the box size
741:         DOUBLE PRECISION :: XJ(3), RIJ(3), MODRIJ, DVDR, D2VDR2, BOX_SIZE(3)688:         DOUBLE PRECISION :: XJ(3), RIJ(3), MODRIJ, DVDR, BOX_SIZE(3)
742: 689: 
743:         ! PER_k: integers for looping over cell possibilities690:         ! PER_k: integers for looping over cell possibilities
744:         ! CELL_RANGE: number of cells in each direction required691:         ! CELL_RANGE: number of cells in each direction required
745:         INTEGER :: PER_X, PER_Y, PER_Z, CELL_RANGE(3)692:         INTEGER :: PER_X, PER_Y, PER_Z, CELL_RANGE(3)
746: 693: 
747: !        DOUBLE PRECISION :: START_ENERGY694: !        DOUBLE PRECISION :: START_ENERGY
748: !        START_ENERGY = ENERGY695: !        START_ENERGY = ENERGY
749: 696: 
750:         BOX_SIZE(1) = BOXLX697:         BOX_SIZE(1) = BOXLX
751:         BOX_SIZE(2) = BOXLY698:         BOX_SIZE(2) = BOXLY
772:                     XJ(3) = COORDS_J(3) + PER_Z719:                     XJ(3) = COORDS_J(3) + PER_Z
773: 720: 
774:                     ! Find the distance and scale by the box size721:                     ! Find the distance and scale by the box size
775:                     RIJ(:) = (COORDS_I(:) - XJ(:))*BOX_SIZE(:)722:                     RIJ(:) = (COORDS_I(:) - XJ(:))*BOX_SIZE(:)
776:                     MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))723:                     MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))
777: 724: 
778:                     IF (MODRIJ .LT. OPP_RCUT) THEN725:                     IF (MODRIJ .LT. OPP_RCUT) THEN
779:                         ! Add energy and gradients726:                         ! Add energy and gradients
780:                         ! Divide by 2 for images of the same atom, to avoid727:                         ! Divide by 2 for images of the same atom, to avoid
781:                         ! double counting728:                         ! double counting
782:                         ! Don't double count energy729:                         ENERGY = ENERGY + MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * &
783:                         IF (IDJ .GE. IDI) THEN730:                                           CALC_ENERGY(MODRIJ)
784:                         ENERGY = ENERGY + & 
785:                             MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * & 
786:                                   CALC_ENERGY(MODRIJ) 
787:                         END IF 
788: 731: 
789:                         IF (GTEST) THEN732:                         IF (GTEST) THEN
790:                             ! Divide gradient by 2 for images of the same atom,733:                             ! Divide gradient by 2 for images of the same atom,
791:                             ! to avoid double counting734:                             ! to avoid double counting
792:                             DVDR = MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * &735:                             DVDR = MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * &
793:                                    CALC_GRAD(MODRIJ) / MODRIJ736:                                    CALC_GRAD(MODRIJ)
794:                             ! We must undo the box size coordinate scaling737:                             ! We must undo the box size coordinate scaling
795:                             GRAD_I(:) = GRAD_I(:) + DVDR * RIJ(:) * BOX_SIZE(:)738:                             GRAD_I(:) = GRAD_I(:) + &
 739:                                         DVDR*RIJ(:)*BOX_SIZE(:)/MODRIJ
 740:                             GRAD_J(:) = GRAD_J(:) - &
 741:                                         DVDR*RIJ(:)*BOX_SIZE(:)/MODRIJ
796:                         END IF ! GTEST742:                         END IF ! GTEST
797:  
798:                         IF (STEST) THEN 
799:                             ! Hessian terms 
800:                             DVDR = MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * & 
801:                                 CALC_GRAD(MODRIJ) / MODRIJ 
802:                             D2VDR2 = MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * & 
803:                                 (CALC_SEC(MODRIJ) - DVDR) / MODRIJ**2 
804:                             ! Same atom, same coordinate 
805:                             HESS(3*IDI-2,3*IDI-2) = HESS(3*IDI-2,3*IDI-2) + & 
806:                                 (D2VDR2 * RIJ(1)**2 + DVDR) * BOX_SIZE(1)**2 
807:                             HESS(3*IDI-1,3*IDI-1) = HESS(3*IDI-1,3*IDI-1) + & 
808:                                 (D2VDR2 * RIJ(2)**2 + DVDR) * BOX_SIZE(2)**2 
809:                             HESS(3*IDI  ,3*IDI  ) = HESS(3*IDI  ,3*IDI  ) + & 
810:                                 (D2VDR2 * RIJ(3)**2 + DVDR) * BOX_SIZE(3)**2 
811:                             ! Same atom, different coordinate 
812:                             HESS(3*IDI-2,3*IDI-1) = HESS(3*IDI-2,3*IDI-1) + & 
813:                                 D2VDR2 * RIJ(1) * RIJ(2) * & 
814:                                 BOX_SIZE(1) * BOX_SIZE(2) 
815:                             HESS(3*IDI-1,3*IDI  ) = HESS(3*IDI-1,3*IDI  ) + & 
816:                                 D2VDR2 * RIJ(2) * RIJ(3) * & 
817:                                 BOX_SIZE(2) * BOX_SIZE(3) 
818:                             HESS(3*IDI  ,3*IDI-2) = HESS(3*IDI  ,3*IDI-2) + & 
819:                                 D2VDR2 * RIJ(3) * RIJ(1) * & 
820:                                 BOX_SIZE(3) * BOX_SIZE(1) 
821:                             HESS(3*IDI-1,3*IDI-2) = HESS(3*IDI-1,3*IDI-2) + & 
822:                                 D2VDR2 * RIJ(2) * RIJ(1) * & 
823:                                 BOX_SIZE(2) * BOX_SIZE(1) 
824:                             HESS(3*IDI  ,3*IDI-1) = HESS(3*IDI  ,3*IDI-1) + & 
825:                                 D2VDR2 * RIJ(3) * RIJ(2) * & 
826:                                 BOX_SIZE(3) * BOX_SIZE(2) 
827:                             HESS(3*IDI-2,3*IDI  ) = HESS(3*IDI-2,3*IDI) + & 
828:                                 D2VDR2 * RIJ(1) * RIJ(3) * & 
829:                                 BOX_SIZE(1) * BOX_SIZE(3) 
830:                             ! Different atom, same coordinate 
831:                             HESS(3*IDI-2,3*IDJ-2) = HESS(3*IDI-2,3*IDJ-2) + & 
832:                                 (- D2VDR2 * RIJ(1)**2 - DVDR) * BOX_SIZE(1)**2 
833:                             HESS(3*IDI-1,3*IDJ-1) = HESS(3*IDI-1,3*IDJ-1) + & 
834:                                 (- D2VDR2 * RIJ(2)**2 - DVDR) * BOX_SIZE(2)**2 
835:                             HESS(3*IDI  ,3*IDJ  ) = HESS(3*IDI  ,3*IDJ  ) + & 
836:                                 (- D2VDR2 * RIJ(3)**2 - DVDR) * BOX_SIZE(3)**2 
837:                             ! Different atom, different coordinate 
838:                             HESS(3*IDI-2,3*IDJ-1) = HESS(3*IDI-2,3*IDJ-1) - & 
839:                                 D2VDR2*RIJ(1)*RIJ(2)*BOX_SIZE(1)*BOX_SIZE(2) 
840:                             HESS(3*IDI-1,3*IDJ  ) = HESS(3*IDI-1,3*IDJ  ) - &  
841:                                 D2VDR2*RIJ(2)*RIJ(3)*BOX_SIZE(2)*BOX_SIZE(3) 
842:                             HESS(3*IDI-2,3*IDJ  ) = HESS(3*IDI-2,3*IDJ  ) - &  
843:                                 D2VDR2*RIJ(3)*RIJ(1)*BOX_SIZE(3)*BOX_SIZE(1) 
844:                             HESS(3*IDJ-1,3*IDI-2) = HESS(3*IDJ-1,3*IDI-2) - & 
845:                                 D2VDR2*RIJ(1)*RIJ(2)*BOX_SIZE(2)*BOX_SIZE(1) 
846:                             HESS(3*IDJ  ,3*IDI-1) = HESS(3*IDJ  ,3*IDI-1) - & 
847:                                 D2VDR2*RIJ(2)*RIJ(3)*BOX_SIZE(3)*BOX_SIZE(2) 
848:                             HESS(3*IDJ  ,3*IDI-2) = HESS(3*IDJ  ,3*IDI-2) - & 
849:                                 D2VDR2*RIJ(3)*RIJ(1)*BOX_SIZE(1)*BOX_SIZE(3) 
850:                         END IF ! End second derivative 
851:                     END IF ! End if less than cutoff743:                     END IF ! End if less than cutoff
852:                 END DO ! z loop744:                 END DO ! z loop
853:             END DO ! y loop745:             END DO ! y loop
854:         END DO ! x loop746:         END DO ! x loop
855: 747: 
856:     END SUBROUTINE OPP_PER748:     END SUBROUTINE OPP_PER
857: 749: 
858: !-------------------------------------------------------------------------------750: !-------------------------------------------------------------------------------
859: !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: !
860: ! Checks whether the unit cell angles are allowed. It is possible to generate832: ! Checks whether the unit cell angles are allowed. It is possible to generate
861: ! sets of angles which do not correspond to a physically possible unit cell.833: ! sets of angles which do not correspond to a physically possible unit cell.
862: ! This would manifest by an imaginary unit cell volume. See834: ! This would manifest by an imaginary unit cell volume. See
863: ! https://journals.iucr.org/a/issues/2011/01/00/au5114/au5114.pdf835: ! https://journals.iucr.org/a/issues/2011/01/00/au5114/au5114.pdf
864: ! 'On the allowed values for the triclinic unit-cell angles',836: ! 'On the allowed values for the triclinic unit-cell angles',
865: ! J. Foadi and G. Evans, Acta Cryst. (2011) A67, 93–95837: ! J. Foadi and G. Evans, Acta Cryst. (2011) A67, 93–95
866: ! ANGLES: the three unit cell angles838: ! ANGLES: the three unit cell angles
867: !839: !
868: !-------------------------------------------------------------------------------840: !-------------------------------------------------------------------------------
869: 841: 
893: ! Finds the energy and gradients for the periodic system with a triclinic865: ! Finds the energy and gradients for the periodic system with a triclinic
894: ! varying unit cell. Optionally with varying potential parameters866: ! varying unit cell. Optionally with varying potential parameters
895: ! IDI, IDJ: id numbers of the particles867: ! IDI, IDJ: id numbers of the particles
896: ! COORDS: coordinates of particles (in crystallographic space)868: ! COORDS: coordinates of particles (in crystallographic space)
897: !         and lattice parameters869: !         and lattice parameters
898: ! GRAD_I: gradients for particle I870: ! GRAD_I: gradients for particle I
899: ! GRAD_J: gradients for particle J871: ! GRAD_J: gradients for particle J
900: ! ENERGY: energy contribution for this pair872: ! ENERGY: energy contribution for this pair
901: ! RIJ: mininum length vector between the particles873: ! RIJ: mininum length vector between the particles
902: ! GTEST: whether to calculate gradients874: ! GTEST: whether to calculate gradients
903: ! STEST: whether to calculate Hessian 
904: ! GRAD_BOX: gradients for the box (angles, then lengths)875: ! GRAD_BOX: gradients for the box (angles, then lengths)
905: ! BD: factors for the lattice derivatives876: ! BD: factors for the lattice derivatives
906: ! GRAD_PARAMS: optional, derivatives of potential parameters877: ! GRAD_PARAMS: optional, derivatives of potential parameters
907: !878: !
908: !-------------------------------------------------------------------------------879: !-------------------------------------------------------------------------------
909: 880: 
910:     SUBROUTINE OPP_PER_TRI(IDI, IDJ, COORDS, GRAD_I, GRAD_J, ENERGY, GTEST, &881:     SUBROUTINE OPP_PER_TRI(IDI, IDJ, COORDS, GRAD_I, GRAD_J, ENERGY, GTEST, &
911:                            STEST, GRAD_BOX, BD, GRAD_PARAMS)882:                            GRAD_BOX, BD, GRAD_PARAMS)
912: 883: 
913:         ! NATOMS: number of coordinates (including lattice parameters)884:         ! NATOMS: number of coordinates (including lattice parameters)
914:         USE COMMONS, ONLY: NATOMS885:         USE COMMONS, ONLY: NATOMS
915: 886: 
916:         ! BOX_DERIV: type for storing box derivative information 
917:         USE TRICLINIC_MOD, ONLY: BOX_DERIV 
918:  
919:         IMPLICIT NONE887:         IMPLICIT NONE
920: 888: 
921:         INTEGER, INTENT(IN) :: IDI, IDJ889:         INTEGER, INTENT(IN) :: IDI, IDJ
922:         DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS)890:         DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS)
923:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_I(3), GRAD_J(3), ENERGY891:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_I(3), GRAD_J(3), ENERGY
924:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_BOX(6)892:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_BOX(6)
925:         LOGICAL, INTENT(IN) :: GTEST, STEST893:         LOGICAL, INTENT(IN) :: GTEST
926:         TYPE(BOX_DERIV), INTENT(IN) :: BD894:         TYPE(BOX_DERIV), INTENT(IN) :: BD
927:         DOUBLE PRECISION, INTENT(INOUT), OPTIONAL :: GRAD_PARAMS(3)895:         DOUBLE PRECISION, INTENT(INOUT), OPTIONAL :: GRAD_PARAMS(3)
928: 896: 
929:         ! RJ: particle coordinates in crystallographic space897:         ! RJ: particle coordinates in crystallographic space
930:         ! RIJ: particle-particle vector in crystallographic space898:         ! RIJ: particle-particle vector in crystallographic space
931:         ! YIJ: particle-particle vector in orthogonal space899:         ! YIJ: particle-particle vector in orthogonal space
932:         ! MODYIJ: particle-particle distance in orthogonal space900:         ! MODYIJ: particle-particle distance in orthogonal space
933:         ! DVDY: derivative of pair potential wrt MODYIJ / MODYIJ901:         ! DVDY: derivative of pair potential wrt MODYIJ
934:         ! TEMP_GRAD: temporary gradient store902:         ! TEMP_GRAD: temporary gradient store
935:         ! TEMP: temporary calculation store 
936:         DOUBLE PRECISION :: RJ(3), RIJ(3), YIJ(3), MODYIJ903:         DOUBLE PRECISION :: RJ(3), RIJ(3), YIJ(3), MODYIJ
937:         DOUBLE PRECISION :: DVDY, TEMP_GRAD(3), TEMP904:         DOUBLE PRECISION :: DVDY, TEMP_GRAD(6)
938: 905: 
939:         ! PER_k: integers for looping over cell possibilities906:         ! PER_k: integers for looping over cell possibilities
940:         ! CELL_RANGE: number of cells in each direction required907:         ! CELL_RANGE: number of cells in each direction required
941:         ! I: loop index908:         ! I, J: loop indices
942:         INTEGER :: PER_X, PER_Y, PER_Z909:         INTEGER :: PER_X, PER_Y, PER_Z
943:         INTEGER :: I910:         INTEGER :: I, J
944: 911: 
945:         ! Loop over different coordinate possibilities.912:         ! Loop over different coordinate possibilities.
946:         DO PER_X = - BD%CELL_RANGE(1), BD%CELL_RANGE(1)913:         DO PER_X = - BD%CELL_RANGE(1), BD%CELL_RANGE(1)
947:             DO PER_Y = - BD%CELL_RANGE(2), BD%CELL_RANGE(2)914:             DO PER_Y = - BD%CELL_RANGE(2), BD%CELL_RANGE(2)
948:                 DO PER_Z = - BD%CELL_RANGE(3), BD%CELL_RANGE(3)915:                 DO PER_Z = - BD%CELL_RANGE(3), BD%CELL_RANGE(3)
949: 916: 
950:                     ! Skip the interaction of a particle with itself in the same917:                     ! Skip the interaction of a particle with itself in the same
951:                     ! cell and don't double count interactions within the unit918:                     ! cell and don't double count interactions within the unit
952:                     ! cell.919:                     ! cell.
953:                     IF (PER_X .EQ. 0 .AND. PER_Y .EQ. 0 .AND. PER_Z .EQ. 0 &920:                     IF (PER_X .EQ. 0 .AND. PER_Y .EQ. 0 .AND. PER_Z .EQ. 0 &
963: 930: 
964:                     ! Find the distance in orthogonal space931:                     ! Find the distance in orthogonal space
965:                     YIJ(:) = MATMUL(BD%ORTHOG(:,:), RIJ(:))932:                     YIJ(:) = MATMUL(BD%ORTHOG(:,:), RIJ(:))
966:                     MODYIJ = SQRT(DOT_PRODUCT(YIJ(:), YIJ(:)))933:                     MODYIJ = SQRT(DOT_PRODUCT(YIJ(:), YIJ(:)))
967: 934: 
968:                     IF (MODYIJ .LT. OPP_RCUT) THEN935:                     IF (MODYIJ .LT. OPP_RCUT) THEN
969:                         ! Add energy and gradients936:                         ! Add energy and gradients
970:                         ! Divide by 2 for images of the same atom, to avoid937:                         ! Divide by 2 for images of the same atom, to avoid
971:                         ! double counting938:                         ! double counting
972:                         ENERGY = ENERGY + MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * &939:                         ENERGY = ENERGY + MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * &
973:                                   CALC_ENERGY(MODYIJ)940:                                           CALC_ENERGY(MODYIJ)
974: 941: 
975:                         ! Gradients 
976:                         IF (GTEST) THEN942:                         IF (GTEST) THEN
977:                             ! Divide gradient by 2 for images of the same atom,943:                             ! Divide gradient by 2 for images of the same atom,
978:                             ! to avoid double counting944:                             ! to avoid double counting
979:                             DVDY = MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * &945:                             DVDY = MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * &
980:                                    CALC_GRAD(MODYIJ) / MODYIJ946:                                    CALC_GRAD(MODYIJ)
981: 947: 
982:                             ! Atom position derivatives948:                             ! Atom position derivatives
983:                             TEMP_GRAD(1:3) = DVDY * &949:                             DO I = 1, 3
984:                                 MATMUL(TRANSPOSE(BD%ORTHOG(:,:)), YIJ(:))950:                                 TEMP_GRAD(I) = DVDY / MODYIJ * &
 951:                                              DOT_PRODUCT(YIJ(:), BD%ORTHOG(:,I))
 952:                             END DO
985:                             GRAD_I(:) = GRAD_I(:) + TEMP_GRAD(1:3)953:                             GRAD_I(:) = GRAD_I(:) + TEMP_GRAD(1:3)
986:                             GRAD_J(:) = GRAD_J(:) - TEMP_GRAD(1:3)954:                             GRAD_J(:) = GRAD_J(:) - TEMP_GRAD(1:3)
987: 955: 
988:                             ! Lattice parameter derivatives956:                             ! Lattice parameter derivatives
 957:                             TEMP_GRAD(:) = 0.D0
989:                             DO I = 1, 6 ! Loop over the different box parameters958:                             DO I = 1, 6 ! Loop over the different box parameters
990:                                 TEMP = DVDY * &959:                                 DO J = 1, 3 ! Loop over x, y, z contributions
991:                                     DOT_PRODUCT(YIJ(:), &960:                                     TEMP_GRAD(I) = TEMP_GRAD(I) + &
992:                                     MATMUL(BD%DERIV(:,:,I), RIJ(:)))961:                                     YIJ(J)*DOT_PRODUCT(BD%DERIV(J,:,I), RIJ(:))
993:                                 GRAD_BOX(I) = GRAD_BOX(I) + TEMP962:                                 END DO
994:                                 IF (I .EQ. 1 .AND. .NOT. STEST) THEN 
995: !                                    WRITE (*, *) 'IDI = ', IDI, 'IDJ = ', IDJ, 'Gradient Contribution = ', TEMP 
996:                                END IF 
997:                             END DO963:                             END DO
 964:                             TEMP_GRAD(:) = TEMP_GRAD(:) * DVDY / MODYIJ
 965:                             GRAD_BOX(:) = GRAD_BOX(:) + TEMP_GRAD(:)
998: 966: 
999:                             ! If necessary, calculate the parameter derivatives967:                             ! If necessary, calculate the parameter derivatives
1000:                             ! Divide gradient by 2 for images of the same atom,968:                             ! Divide gradient by 2 for images of the same atom,
1001:                             ! to avoid double counting969:                             ! to avoid double counting
1002:                             IF(PRESENT(GRAD_PARAMS) .AND. OPP_MODE .EQ. 2)&970:                             IF(PRESENT(GRAD_PARAMS) .AND. OPP_MODE .EQ. 2)&
1003:                                 GRAD_PARAMS(:) = GRAD_PARAMS(:) + &971:                                 GRAD_PARAMS(:) = GRAD_PARAMS(:) + &
1004:                                     MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * &972:                                     MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * &
1005:                                     CALC_GRAD_PARAMS(MODYIJ)973:                                     CALC_GRAD_PARAMS(MODYIJ)
1006: 974: 
1007:                         END IF ! GTEST975:                         END IF ! GTEST
1008:  
1009:                         IF (STEST) THEN 
1010:                             CALL HESS_CONTRIB(RIJ(:), BD, IDI, IDJ) 
1011:                         END IF ! End second derivative 
1012:                     END IF ! End if less than cutoff976:                     END IF ! End if less than cutoff
1013:                 END DO ! z loop977:                 END DO ! z loop
1014:             END DO ! y loop978:             END DO ! y loop
1015:         END DO ! x loop979:         END DO ! x loop
1016: 980: 
1017:     END SUBROUTINE OPP_PER_TRI981:     END SUBROUTINE OPP_PER_TRI
1018: 982: 
1019: !-------------------------------------------------------------------------------983: !-------------------------------------------------------------------------------
1020: !984: !
1021: ! Calculates the second derivative contributions 
1022: ! RIJ: vector between particles in fractional coordinates 
1023: ! BD: box derivative information 
1024: ! IDI, IDJ: indices of the atom pair 
1025: ! 
1026: !------------------------------------------------------------------------------- 
1027:     SUBROUTINE HESS_CONTRIB(RIJ, BD, IDI, IDJ) 
1028:  
1029:         ! Number of atoms, including unit cell parameters 
1030:         USE COMMONS, ONLY: NATOMS 
1031:  
1032:         ! HESS: the Hessian matrix 
1033:         USE MODHESS, ONLY: HESS 
1034:  
1035:         ! BOX_DERIV: type for storing box derivative information 
1036:         USE TRICLINIC_MOD, ONLY: BOX_DERIV 
1037:  
1038:         IMPLICIT NONE 
1039:  
1040:         DOUBLE PRECISION, INTENT(IN) :: RIJ(3) 
1041:         TYPE(BOX_DERIV), INTENT(IN) :: BD 
1042:         INTEGER, INTENT(IN) :: IDI, IDJ 
1043:  
1044:         ! YIJ: particle-particle vector in Cartesian space 
1045:         ! MODYIJ: particle-particle distance in Cartesian space 
1046:         ! DVDY: derivative of pair potential wrt MODYIJ 
1047:         ! D2VDY2: D(DVDY/MODYIJ)/DMODYIJ * 1/MODYIJ 
1048:         ! TEMP, TEMP_GRAD: temporary calculation store 
1049:         DOUBLE PRECISION :: YIJ(3), MODYIJ, DVDY, D2VDY2, TEMP, TEMP_GRAD(3) 
1050:  
1051:         ! I, J: loop indices 
1052:         INTEGER :: I, J 
1053:  
1054:         ! Calculate some factors 
1055:         YIJ(:) = MATMUL(BD%ORTHOG(:,:), RIJ(:)) 
1056:         MODYIJ = SQRT(DOT_PRODUCT(YIJ(:), YIJ(:))) 
1057:         DVDY = MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * CALC_GRAD(MODYIJ) / MODYIJ 
1058:         D2VDY2 = (MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * CALC_SEC(MODYIJ) - DVDY) & 
1059:             / MODYIJ**2 
1060:  
1061:         ! Pure atom translations 
1062:         ! Loop over the atom coordinates twice 
1063:         DO I = 1, 3 
1064:             DO J = 1, 3 
1065:                 TEMP = D2VDY2 * DOT_PRODUCT(YIJ(:),BD%ORTHOG(:,I)) * & 
1066:                     DOT_PRODUCT(YIJ(:),BD%ORTHOG(:,J)) + DVDY * BD%ORTHOGSQ(I,J) 
1067:                 ! Same atom, same or different coord 
1068:                 HESS(3*IDI-3+I,3*IDI-3+J) = HESS(3*IDI-3+I,3*IDI-3+J) + TEMP 
1069:                 HESS(3*IDJ-3+I,3*IDJ-3+J) = HESS(3*IDJ-3+I,3*IDJ-3+J) + TEMP 
1070:                 ! Different atom, same or different coord 
1071:                 HESS(3*IDI-3+I,3*IDJ-3+J) = HESS(3*IDI-3+I,3*IDJ-3+J) - TEMP 
1072:                 HESS(3*IDJ-3+J,3*IDI-3+I) = HESS(3*IDJ-3+J,3*IDI-3+I) - TEMP 
1073:             END DO 
1074:         END DO 
1075:  
1076:         ! Mixed atom and unit cell parameter 
1077:         DO I = 1, 6 ! loop over unit cell parameters 
1078:             TEMP_GRAD(:) = D2VDY2 * & 
1079:                 DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,I), RIJ(:))) * & 
1080:                 MATMUL(TRANSPOSE(BD%ORTHOG(:,:)), YIJ(:)) + & 
1081:                 DVDY * (MATMUL(TRANSPOSE(BD%DERIV(:,:,I)), YIJ(:)) + & 
1082:                 MATMUL(TRANSPOSE(BD%ORTHOG(:,:)), & 
1083:                 MATMUL(BD%DERIV(:,:,I), RIJ(:)))) 
1084:             HESS(3*IDI-2:3*IDI,3*NATOMS-6+I) = & 
1085:                 HESS(3*IDI-2:3*IDI,3*NATOMS-6+I) + TEMP_GRAD(:) 
1086:             HESS(3*NATOMS-6+I,3*IDI-2:3*IDI) = & 
1087:                 HESS(3*NATOMS-6+I,3*IDI-2:3*IDI) + TEMP_GRAD(:) 
1088:             HESS(3*IDJ-2:3*IDJ,3*NATOMS-6+I) = & 
1089:                 HESS(3*IDJ-2:3*IDJ,3*NATOMS-6+I) - TEMP_GRAD(:) 
1090:             HESS(3*NATOMS-6+I,3*IDJ-2:3*IDJ) = & 
1091:                 HESS(3*NATOMS-6+I,3*IDJ-2:3*IDJ) - TEMP_GRAD(:) 
1092:         END DO ! end loop over unit cell parameters 
1093:  
1094:         ! Pure unit cell parameters 
1095:         DO I = 1, 6 ! outer loop over cell parameters 
1096:             ! Diagonals 
1097:             TEMP = D2VDY2 * & 
1098:                 DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,I), RIJ(:)))**2 + & 
1099:                 DVDY * (DOT_PRODUCT(MATMUL(BD%DERIV(:,:,I), RIJ(:)), & 
1100:                 MATMUL(BD%DERIV(:,:,I), RIJ(:))) + & 
1101:                 DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV2(:,:,I,I), RIJ(:)))) 
1102:             HESS(3*NATOMS-6+I,3*NATOMS-6+I) = HESS(3*NATOMS-6+I,3*NATOMS-6+I) + TEMP 
1103:             IF (I .EQ. 1) THEN 
1104: !                WRITE (*, *) 'IDI = ', IDI, 'IDJ = ', IDJ, 'Hessian contribution = ', TEMP 
1105: !                WRITE (*, *) 'MODYIJ = ', MODYIJ 
1106: !                WRITE (*, *) 'RIJ = ', RIJ 
1107: !                WRITE (*, *) 'YIJ = ', YIJ 
1108: !                WRITE (*, *) 'CALC_GRAD(MODYIJ) = ', CALC_GRAD(MODYIJ) 
1109: !                WRITE (*, *) 'DVDY =', DVDY 
1110: !                WRITE (*, *) 'D2VDY2 = ', D2VDY2 
1111:             END IF 
1112:             ! Off diagonals 
1113:             DO J = I + 1, 6 ! inner loop over cell parameters 
1114:                 TEMP = D2VDY2 * & 
1115:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,I), RIJ(:))) * & 
1116:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,J), RIJ(:))) + & 
1117:                     DVDY * (DOT_PRODUCT(MATMUL(BD%DERIV(:,:,I), RIJ(:)), & 
1118:                     MATMUL(BD%DERIV(:,:,J), RIJ(:))) + & 
1119:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV2(:,:,I,J), RIJ(:)))) 
1120:                 HESS(3*NATOMS-6+I,3*NATOMS-6+J) = & 
1121:                     HESS(3*NATOMS-6+I,3*NATOMS-6+J) + TEMP 
1122:                 HESS(3*NATOMS-6+J,3*NATOMS-6+I) = & 
1123:                     HESS(3*NATOMS-6+J,3*NATOMS-6+I) + TEMP 
1124:             END DO ! end inner loop over cell parameters 
1125:         END DO ! end outer loop over cell parameters 
1126:  
1127:         IF (OPP_MODE .EQ. 2) THEN 
1128:             ! Potential parameter contributions 
1129:             CALL CALC_HESS_PARAMS(RIJ, BD, IDI, IDJ) 
1130:         END IF 
1131:  
1132:     END SUBROUTINE HESS_CONTRIB 
1133:  
1134: !------------------------------------------------------------------------------- 
1135: ! 
1136: ! Calculates the contribution to the derivatives of the potential wrt the985: ! Calculates the contribution to the derivatives of the potential wrt the
1137: ! potential parameters OPP_K and OPP_PHI986: ! potential parameters OPP_K and OPP_PHI
1138: ! MODYIJ: distance between the particles for this contribution987: ! MODYIJ: distance between the particles for this contribution
1139: !988: !
1140: !-------------------------------------------------------------------------------989: !-------------------------------------------------------------------------------
1141: 990: 
1142:     PURE FUNCTION CALC_GRAD_PARAMS(MODYIJ)991:     PURE FUNCTION CALC_GRAD_PARAMS(MODYIJ)
1143: 992: 
1144:         IMPLICIT NONE993:         IMPLICIT NONE
1145: 994: 
1159:         C = COS(OPP_K * (OPP_RCUT - 1) - OPP_PHI)1008:         C = COS(OPP_K * (OPP_RCUT - 1) - OPP_PHI)
1160: 1009: 
1161:         IF (BTEST(OPP_PARAMS, 0)) THEN1010:         IF (BTEST(OPP_PARAMS, 0)) THEN
1162:             ! OPP_K derivative1011:             ! OPP_K derivative
1163:             ! Unmodified potential term1012:             ! Unmodified potential term
1164:             T(1) = - SIN(OPP_K * (MODYIJ-1) - OPP_PHI) * (MODYIJ-1) / MODYIJ**31013:             T(1) = - SIN(OPP_K * (MODYIJ-1) - OPP_PHI) * (MODYIJ-1) / MODYIJ**3
1165:             ! Term to make the potential go to zero at OPP_RCUT1014:             ! Term to make the potential go to zero at OPP_RCUT
1166:             T(1) = T(1) + S * (OPP_RCUT - 1) / OPP_RCUT **31015:             T(1) = T(1) + S * (OPP_RCUT - 1) / OPP_RCUT **3
1167:             ! Terms to make the first derivative go to zero at OPP_RCUT1016:             ! Terms to make the first derivative go to zero at OPP_RCUT
1168:             T(1) = T(1) + (MODYIJ - OPP_RCUT) * &1017:             T(1) = T(1) + (MODYIJ - OPP_RCUT) * &
1169:                 (C * OPP_K * (OPP_RCUT - 1) + S - &1018:                 (C * OPP_K * (OPP_RCUT - 1) + S  - &
1170:                  S * 3.D0 * (OPP_RCUT - 1) / OPP_RCUT ) / OPP_RCUT**31019:                  S * 3.D0 * (OPP_RCUT - 1) / OPP_RCUT ) / OPP_RCUT**3
1171:             ! Terms to make the second derivative go to zero at OPP_RCUT1020:             ! Terms to make the second derivative go to zero at OPP_RCUT
1172: !            T(1) = T(1) - 0.5D0 * (MODYIJ - OPP_RCUT)**2 * &1021: !            T(1) = T(1) - 0.5D0 * (MODYIJ - OPP_RCUT)**2 * &
1173: !                   (S * OPP_K**2 * (OPP_RCUT - 1) / OPP_RCUT**3 - &1022: !                   (S * OPP_K**2 * (OPP_RCUT - 1) / OPP_RCUT**3 - &
1174: !                    C * 2.D0 * OPP_K / OPP_RCUT**3 + &1023: !                    C * 2.D0 * OPP_K / OPP_RCUT**3 + &
1175: !                    C * 6.D0 * OPP_K * (OPP_RCUT - 1) / OPP_RCUT**4 + &1024: !                    C * 6.D0 * OPP_K * (OPP_RCUT - 1) / OPP_RCUT**4 + &
1176: !                    S * 6.D0 / OPP_RCUT**4 - &1025: !                    S * 6.D0 / OPP_RCUT**4 - &
1177: !                    S * 12.D0 * (OPP_RCUT - 1) / OPP_RCUT**5)1026: !                    S * 12.D0 * (OPP_RCUT - 1) / OPP_RCUT**5)
1178:         END IF1027:         END IF
1179:         IF (BTEST(OPP_PARAMS, 1)) THEN1028:         IF (BTEST(OPP_PARAMS, 1)) THEN
1190: !                   ( - S * OPP_K**2 / OPP_RCUT**3 - &1039: !                   ( - S * OPP_K**2 / OPP_RCUT**3 - &
1191: !                    C * 6.D0 * OPP_K / OPP_RCUT**4 + S * 12.D0 / OPP_RCUT**5)           1040: !                    C * 6.D0 * OPP_K / OPP_RCUT**4 + S * 12.D0 / OPP_RCUT**5)           
1192:         END IF1041:         END IF
1193: 1042: 
1194:         CALC_GRAD_PARAMS(:) = T(:)1043:         CALC_GRAD_PARAMS(:) = T(:)
1195: 1044: 
1196:     END FUNCTION CALC_GRAD_PARAMS1045:     END FUNCTION CALC_GRAD_PARAMS
1197: 1046: 
1198: !-------------------------------------------------------------------------------1047: !-------------------------------------------------------------------------------
1199: !1048: !
1200: ! Calculates the second derivative contributions for the potential parameters 
1201: ! RIJ: vector between particles in fractional coordinates 
1202: ! BD: box derivative information 
1203: ! IDI, IDJ: indices of the atom pair 
1204: ! 
1205: !------------------------------------------------------------------------------- 
1206:  
1207:     SUBROUTINE CALC_HESS_PARAMS(RIJ, BD, IDI, IDJ) 
1208:  
1209:         ! Number of atoms, including unit cell parameters 
1210:         USE COMMONS, ONLY: NATOMS 
1211:  
1212:         ! HESS: the Hessian matrix 
1213:         USE MODHESS, ONLY: HESS 
1214:  
1215:         ! BOX_DERIV: type for storing box derivative information 
1216:         USE TRICLINIC_MOD, ONLY: BOX_DERIV 
1217:  
1218:         IMPLICIT NONE 
1219:  
1220:         DOUBLE PRECISION, INTENT(IN) :: RIJ(3) 
1221:         TYPE(BOX_DERIV), INTENT(IN) :: BD 
1222:         INTEGER, INTENT(IN) :: IDI, IDJ 
1223:  
1224:         ! YIJ: particle-particle vector in Cartesian space 
1225:         ! MODYIJ: particle-particle distance in Cartesian space 
1226:         ! S, C, SC, CC: recurring sin and cos terms 
1227:         ! TEMP, TEMP_GRAD: temporary calculation stores 
1228:         ! D2VDK2, D2VDPHI2, D2VDKDPHI: pure parameter second derivatives 
1229:         ! D2VDYDK, D2VDYDPHI; mixed second derivatives of pair potential 
1230:         DOUBLE PRECISION :: YIJ(3), MODYIJ, S, C, SC, CC, TEMP, TEMP_GRAD(3) 
1231:         DOUBLE PRECISION :: D2VDK2, D2VDPHI2, D2VDKDPHI, D2VDYDK, D2VDYDPHI 
1232:  
1233:         ! I: loop index 
1234:         INTEGER :: I 
1235:  
1236:         ! Calculate some factors 
1237:         YIJ(:) = MATMUL(BD%ORTHOG(:,:), RIJ(:)) 
1238:         MODYIJ = SQRT(DOT_PRODUCT(YIJ(:), YIJ(:))) 
1239:         S = SIN(OPP_K * (MODYIJ - 1) - OPP_PHI) 
1240:         C = COS(OPP_K * (MODYIJ - 1) - OPP_PHI) 
1241:         SC = SIN(OPP_K * (OPP_RCUT - 1) - OPP_PHI) 
1242:         CC = COS(OPP_K * (OPP_RCUT - 1) - OPP_PHI) 
1243:  
1244:         IF (BTEST(OPP_PARAMS, 0)) THEN 
1245:             ! OPP_K pure derivative 
1246:             ! Unmodified potential term 
1247:             D2VDK2 = - C * (MODYIJ - 1)**2 / MODYIJ**3 
1248:             ! Term to make the potential go to zero at OPP_RCUT 
1249:             D2VDK2 = D2VDK2 + CC * (OPP_RCUT - 1)**2 / OPP_RCUT **3 
1250:             ! Terms to make the first derivative go to zero at OPP_RCUT 
1251:             D2VDK2 = D2VDK2 + (MODYIJ - OPP_RCUT) * & 
1252:                 (- SC * OPP_K * (OPP_RCUT-1)**2 + CC * 2.D0 * (OPP_RCUT-1) - & 
1253:                 CC * 3.D0 * (OPP_RCUT - 1)**2 / OPP_RCUT ) / OPP_RCUT**3 
1254:             D2VDK2 = D2VDK2 * MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) 
1255:             HESS(3*NATOMS-8,3*NATOMS-8) = HESS(3*NATOMS-8,3*NATOMS-8) + D2VDK2 
1256:  
1257:             ! OPP_K and position mixed derivative 
1258:             ! Unmodified potential term 
1259:             D2VDYDK = (- C * OPP_K * (MODYIJ - 1) - S + & 
1260:                 S * 3.D0 * (MODYIJ - 1) / MODYIJ) / MODYIJ**3 
1261:             ! Terms to make the first derivative go to zero at OPP_RCUT 
1262:             D2VDYDK = D2VDYDK + (CC * OPP_K * (OPP_RCUT - 1) + SC - & 
1263:                 SC * 3.D0 * (OPP_RCUT - 1) / OPP_RCUT) / OPP_RCUT**3 
1264:             D2VDYDK = D2VDYDK * MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) / MODYIJ 
1265:             TEMP_GRAD(:) = D2VDYDK * MATMUL(TRANSPOSE(BD%ORTHOG(:,:)), YIJ(:)) 
1266:             HESS(3*IDI-2:3*IDI,3*NATOMS-8) = & 
1267:                 HESS(3*IDI-2:3*IDI,3*NATOMS-8) + TEMP_GRAD(:) 
1268:             HESS(3*NATOMS-8,3*IDI-2:3*IDI) = & 
1269:                 HESS(3*NATOMS-8,3*IDI-2:3*IDI) + TEMP_GRAD(:) 
1270:             HESS(3*IDJ-2:3*IDJ,3*NATOMS-8) = & 
1271:                 HESS(3*IDJ-2:3*IDJ,3*NATOMS-8) - TEMP_GRAD(:) 
1272:             HESS(3*NATOMS-8,3*IDJ-2:3*IDJ) = & 
1273:                 HESS(3*NATOMS-8,3*IDJ-2:3*IDJ) - TEMP_GRAD(:) 
1274:  
1275:             ! OPP_K and unit cell parameter mixed derivative 
1276:             DO I = 1, 6 ! loop over unit cell parameters 
1277:                 TEMP = D2VDYDK * & 
1278:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,I), RIJ(:))) 
1279:                 HESS(3*NATOMS-6+I,3*NATOMS-8) = & 
1280:                     HESS(3*NATOMS-6+I,3*NATOMS-8) + TEMP 
1281:                 HESS(3*NATOMS-8,3*NATOMS-6+I) = & 
1282:                     HESS(3*NATOMS-8,3*NATOMS-6+I) + TEMP 
1283:             END DO ! end loop over unit cell parameters 
1284:         END IF         
1285:  
1286:         IF (BTEST(OPP_PARAMS, 1)) THEN 
1287:             ! OPP_PHI pure derivative 
1288:             ! Unmodified potential term 
1289:             D2VDPHI2 = - C / MODYIJ**3 
1290:             ! Term to make the potential go to zero at OPP_RCUT 
1291:             D2VDPHI2 = D2VDPHI2 + CC / OPP_RCUT **3 
1292:             ! Terms to make the first derivative go to zero at OPP_RCUT 
1293:             D2VDPHI2 = D2VDPHI2 + (MODYIJ - OPP_RCUT) * & 
1294:                 (- SC * OPP_K - CC * 3.D0 / OPP_RCUT ) / OPP_RCUT**3 
1295:             D2VDPHI2 = D2VDPHI2 * MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) 
1296:             HESS(3*NATOMS-7,3*NATOMS-7) = HESS(3*NATOMS-7,3*NATOMS-7) + D2VDPHI2 
1297:  
1298:             ! OPP_PHI and position mixed derivative 
1299:             ! Unmodified potential term 
1300:             D2VDYDPHI = (C * OPP_K - S * 3.D0 / MODYIJ) / MODYIJ**3 
1301:             ! Terms to make the first derivative go to zero at OPP_RCUT 
1302:             D2VDYDPHI = D2VDYDPHI + ( - CC * OPP_K + & 
1303:                 SC * 3.D0 / OPP_RCUT) / OPP_RCUT**3 
1304:             D2VDYDPHI = D2VDYDPHI * MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) / MODYIJ 
1305:             TEMP_GRAD(:) = D2VDYDPHI * MATMUL(TRANSPOSE(BD%ORTHOG(:,:)), YIJ(:)) 
1306:             HESS(3*IDI-2:3*IDI,3*NATOMS-7) = & 
1307:                 HESS(3*IDI-2:3*IDI,3*NATOMS-7) + TEMP_GRAD(:) 
1308:             HESS(3*NATOMS-7,3*IDI-2:3*IDI) = & 
1309:                 HESS(3*NATOMS-7,3*IDI-2:3*IDI) + TEMP_GRAD(:) 
1310:             HESS(3*IDJ-2:3*IDJ,3*NATOMS-7) = & 
1311:                 HESS(3*IDJ-2:3*IDJ,3*NATOMS-7) - TEMP_GRAD(:) 
1312:             HESS(3*NATOMS-7,3*IDJ-2:3*IDJ) = & 
1313:                 HESS(3*NATOMS-7,3*IDJ-2:3*IDJ) - TEMP_GRAD(:) 
1314:  
1315:             ! OPP_K and unit cell parameter mixed derivative 
1316:             DO I = 1, 6 ! loop over unit cell parameters 
1317:                 TEMP = D2VDYDPHI * & 
1318:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,I), RIJ(:))) 
1319:                 HESS(3*NATOMS-6+I,3*NATOMS-7) = & 
1320:                     HESS(3*NATOMS-6+I,3*NATOMS-7) + TEMP 
1321:                 HESS(3*NATOMS-7,3*NATOMS-6+I) = & 
1322:                     HESS(3*NATOMS-7,3*NATOMS-6+I) + TEMP 
1323:             END DO ! end loop over unit cell parameters 
1324:         END IF 
1325:  
1326:         IF (BTEST(OPP_PARAMS, 0) .AND. BTEST(OPP_PARAMS, 1)) THEN 
1327:             ! OPP_K and OPP_PHI derivative 
1328:             ! Unmodified potential term 
1329:             D2VDKDPHI = C * (MODYIJ - 1) / MODYIJ**3 
1330:             ! Term to make the potential go to zero at OPP_RCUT 
1331:             D2VDKDPHI = D2VDKDPHI - CC * (OPP_RCUT - 1) / OPP_RCUT**3 
1332:             ! Terms to make the first derivative go to zero at OPP_RCUT 
1333:             D2VDKDPHI = D2VDKDPHI + (MODYIJ - OPP_RCUT) * & 
1334:                 (SC * OPP_K * (OPP_RCUT - 1) -  CC + & 
1335:                 CC * 3.D0 * (OPP_RCUT - 1) / OPP_RCUT) / OPP_RCUT**3 
1336:             D2VDKDPHI = D2VDKDPHI * MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) 
1337:             HESS(3*NATOMS-8,3*NATOMS-7) = HESS(3*NATOMS-8,3*NATOMS-7)+D2VDKDPHI 
1338:             HESS(3*NATOMS-7,3*NATOMS-8) = HESS(3*NATOMS-7,3*NATOMS-8)+D2VDKDPHI 
1339:         END IF 
1340:  
1341:     END SUBROUTINE CALC_HESS_PARAMS 
1342:  
1343: !------------------------------------------------------------------------------- 
1344: ! 
1345: ! Repel the unit cell volume from approaching 0, which repels the unit cell from1049: ! Repel the unit cell volume from approaching 0, which repels the unit cell from
1346: ! adopting physically impossible angle combinations. Uses a WCA potential1050: ! adopting physically impossible angle combinations. Uses a WCA potential
1347: ! ENERGY: energy of the system1051: ! ENERGY: energy of the system
1348: ! GRAD_ANGLES: gradients for the box angles1052: ! GRAD_ANGLES: gradients for the box angles
1349: ! BD: information on the box and its derivatives1053: ! BD: information on the box and its derivatives
1350: ! GTEST: whether to calculate gradients1054: ! GTEST: whether to calculate gradients
1351: ! STEST: whether to calculate second derivatives 
1352: !1055: !
1353: !-------------------------------------------------------------------------------1056: !-------------------------------------------------------------------------------
1354: 1057: 
1355:     SUBROUTINE CONSTRAIN_VOLUME(ENERGY, GRAD_ANGLES, BD, GTEST, STEST)1058:     SUBROUTINE CONSTRAIN_VOLUME(ENERGY, GRAD_ANGLES, BD, GTEST)
1356:  
1357:         ! NATOMS: number of atoms, including unit cell parameters 
1358:         USE COMMONS, ONLY: NATOMS 
1359:  
1360:         ! HESS: the Hessian matrix 
1361:         USE MODHESS, ONLY: HESS 
1362:  
1363:         ! BOX_DERIV: type for storing box derivative information 
1364:         USE TRICLINIC_MOD, ONLY: BOX_DERIV 
1365: 1059: 
1366:         IMPLICIT NONE1060:         IMPLICIT NONE
1367: 1061: 
1368:         DOUBLE PRECISION, INTENT(INOUT) :: ENERGY, GRAD_ANGLES(3)1062:         DOUBLE PRECISION, INTENT(INOUT) :: ENERGY, GRAD_ANGLES(3)
1369:         TYPE(BOX_DERIV), INTENT(IN) :: BD1063:         TYPE(BOX_DERIV), INTENT(IN) :: BD
1370:         LOGICAL, INTENT(IN) :: GTEST, STEST1064:         LOGICAL, INTENT(IN) :: GTEST
1371:  
1372:         ! DEDV: derivative of energy wrt volume 
1373:         ! D2EDV2: second derivative of energy wrt volume 
1374:         DOUBLE PRECISION :: DEDV, D2EDV2 
1375:  
1376:         ! I: loop index 
1377:         INTEGER :: I 
1378: 1065: 
1379:         ! Add a purely repulsive (away from zero) WCA energy term1066:         ! Add a purely repulsive (away from zero) WCA energy term
1380:         IF (BD%V .LT. V_SIGMA * 2.D0**(1.D0/6.D0)) THEN1067:         IF (BD%V .LT. V_SIGMA**(1.D0/6.D0)) THEN
1381:             ENERGY = ENERGY + 4.D0 * V_EPS * &1068:             ENERGY = ENERGY + 4.D0 * V_EPS * &
1382:                      ((V_SIGMA / BD%V)**(12) - (V_SIGMA / BD%V)**(6)) + V_EPS1069:                      ((V_SIGMA / BD%V)**(12) - (V_SIGMA / BD%V)**(6)) + V_EPS
1383: 1070: 
1384:             IF (GTEST) THEN1071:             IF (GTEST) THEN
1385:                 ! Add the gradients, if necessary1072:                 ! Add the gradients, if necessary
1386:                 DEDV = 24.D0 * V_EPS / V_SIGMA * &1073:                 GRAD_ANGLES(:) = GRAD_ANGLES(:) + 24.D0 * V_EPS / V_SIGMA * &
1387:                     ((V_SIGMA / BD%V)**(7) - 2.D0 * (V_SIGMA / BD%V)**(13))1074:                     ((V_SIGMA / BD%V)**(7) - 2.D0 * (V_SIGMA / BD%V)**(13)) * &
1388:                 GRAD_ANGLES(:) = GRAD_ANGLES(:) +  DEDV * BD%V_DERIV(:)1075:                     BD%V_DERIV(:)
1389:             END IF 
1390:  
1391:             IF (STEST) THEN 
1392:                 ! Add the Hessian contributions, if necessary 
1393:                 DEDV = 24.D0 * V_EPS / V_SIGMA * & 
1394:                     ((V_SIGMA / BD%V)**(7) - 2.D0 * (V_SIGMA / BD%V)**(13)) 
1395:                 D2EDV2 = 24.D0 * V_EPS / V_SIGMA**2 * & 
1396:                     (26.D0 * (V_SIGMA/BD%V)**(14) - 7.D0 * (V_SIGMA/BD%V)**(8)) 
1397:  
1398:                 DO I = 1, 3 ! outer loop over angles 
1399:                     HESS(3*NATOMS-5+I,3*NATOMS-5:3*NATOMS-3) = & 
1400:                         HESS(3*NATOMS-5+I,3*NATOMS-5:3*NATOMS-3) + & 
1401:                         D2EDV2 * BD%V_DERIV(I) * BD%V_DERIV(:) + & 
1402:                         DEDV * BD%V_DERIV2(I, :) 
1403:                 END DO ! end outer loop over angles 
1404:             END IF1076:             END IF
1405:         END IF1077:         END IF
1406: 1078: 
1407:     END SUBROUTINE CONSTRAIN_VOLUME1079:     END SUBROUTINE CONSTRAIN_VOLUME
1408: 1080: 
1409: !-------------------------------------------------------------------------------1081: !-------------------------------------------------------------------------------
1410: !1082: !
1411: ! Impose limits on the allowed range of the potential parameters OPP_K by adding1083: ! Impose limits on the allowed range of the potential parameters OPP_K by adding
1412: ! a harmonic repulsion outside the allowed range.1084: ! a harmonic repulsion outside the allowed range.
1413: ! ENERGY: energy of the system1085: ! ENERGY: energy of the system
1414: ! GRAD_PARAMS: gradients for the potential parameters1086: ! GRAD_PARAMS: gradients for the potential parameters
1415: ! GTEST: whether to calculate gradients1087: ! GTEST: whether to calculate gradients
1416: ! STEST: whether to calculate second derivatives 
1417: !1088: !
1418: !-------------------------------------------------------------------------------1089: !-------------------------------------------------------------------------------
1419:     SUBROUTINE CONSTRAIN_PARAMETERS(ENERGY, GRAD_PARAMS, GTEST, STEST)1090:     SUBROUTINE CONSTRAIN_PARAMETERS(ENERGY, GRAD_PARAMS, GTEST)
1420:  
1421:         ! NATOMS: number of atoms, including unit cell and potential parameters 
1422:         USE COMMONS, ONLY: NATOMS 
1423:  
1424:         ! HESS: the Hessian matrix 
1425:         USE MODHESS, ONLY: HESS         
1426: 1091: 
1427:         IMPLICIT NONE1092:         IMPLICIT NONE
1428: 1093: 
1429:         DOUBLE PRECISION, INTENT(INOUT) :: ENERGY, GRAD_PARAMS(3)1094:         DOUBLE PRECISION, INTENT(INOUT) :: ENERGY, GRAD_PARAMS(3)
1430:         LOGICAL, INTENT(IN) :: GTEST, STEST1095:         LOGICAL, INTENT(IN) :: GTEST
1431: 1096: 
1432:         ! EN_CON: constraint contribution to the energy1097:         ! EN_CON: constraint contribution to the energy
1433:         DOUBLE PRECISION :: EN_CON1098:         DOUBLE PRECISION :: EN_CON
1434: 1099: 
1435:         ! Initialise output variable1100:         ! Initialise output variable
1436:         EN_CON = 0.D01101:         EN_CON = 0.D0
1437: 1102: 
1438:         ! OPP_K energy contribution1103:         ! OPP_K energy contribution
1439:         IF (BTEST(OPP_PARAMS, 0)) THEN1104:         IF (BTEST(OPP_PARAMS, 0)) THEN
1440:             EN_CON = EN_CON + 0.5D0 * K_REPULSION * &1105:             EN_CON = EN_CON + 0.5D0 * K_REPULSION * &
1441:                 (MERGE(OPP_K - K_UPPER, 0.D0, &1106:                 (MERGE(OPP_K - K_UPPER, 0.D0, &
1442:                        OPP_K .GT. K_UPPER)**2 + &1107:                        OPP_K .GT. K_UPPER)**2 + &
1443:                  MERGE(OPP_K - K_LOWER, 0.D0, &1108:                  MERGE(OPP_K - K_LOWER, 0.D0, &
1444:                        OPP_K .LT. K_LOWER)**2)1109:                        OPP_K .LT. K_LOWER)**2)
1445:  
1446:             IF (GTEST) THEN1110:             IF (GTEST) THEN
1447:                 ! OPP_K gradient contribution1111:             ! OPP_K gradient contribution
1448:                 GRAD_PARAMS(1) = GRAD_PARAMS(1) + K_REPULSION * &1112:                 GRAD_PARAMS(1) = GRAD_PARAMS(1) + K_REPULSION * &
1449:                     (MERGE(OPP_K - K_UPPER, 0.D0, &1113:                     (MERGE(OPP_K - K_UPPER, 0.D0, &
1450:                            OPP_K .GT. K_UPPER) + &1114:                            OPP_K .GT. K_UPPER) + &
1451:                      MERGE(OPP_K - K_LOWER, 0.D0, &1115:                      MERGE(OPP_K - K_LOWER, 0.D0, &
1452:                            OPP_K .LT. K_LOWER))1116:                            OPP_K .LT. K_LOWER))
1453:             END IF1117:             END IF
1454:  
1455:             IF (STEST) THEN 
1456:                 ! OPP_K Hessian contribution 
1457:                 HESS(3*NATOMS-8,3*NATOMS-8) = HESS(3*NATOMS-8,3*NATOMS-8) + & 
1458:                     MERGE(K_REPULSION, 0.D0, OPP_K .GT. K_UPPER) + & 
1459:                     MERGE(K_REPULSION, 0.D0, OPP_K .LT. K_LOWER) 
1460:             END IF  
1461:         END IF1118:         END IF
1462: 1119: 
1463:         ! Add the energy contribution1120:         ! Add the energy contribution
1464:         ENERGY = ENERGY + EN_CON1121:         ENERGY = ENERGY + EN_CON
1465: 1122: 
1466:     END SUBROUTINE CONSTRAIN_PARAMETERS1123:     END SUBROUTINE CONSTRAIN_PARAMETERS
1467: 1124: 
1468: !-------------------------------------------------------------------------------1125: !-------------------------------------------------------------------------------
1469: !1126: !
1470: ! Outputs the coordinates to the file opp.xyz1127: ! Outputs the coordinates to the file ljgauss.xyz
1471: !1128: !
1472: !-------------------------------------------------------------------------------1129: !-------------------------------------------------------------------------------
1473:     SUBROUTINE VIEW_OPP()1130:     SUBROUTINE VIEW_OPP()
1474: 1131: 
1475:         ! NATOMS: number of coordinates1132:         ! NATOMS: number of coordinates
1476:         ! NSAVE: number of minima to write1133:         ! NSAVE: number of minima to write
1477:         ! BOXx: Box dimensions1134:         ! BOXx: Box dimensions
1478:         ! PERIODIC: whether periodic boundary conditions are in use1135:         ! PERIODIC: whether periodic boundary conditions are in use
1479:         USE COMMONS, ONLY: NATOMS, NSAVE, PERIODIC, BOXLX, BOXLY, BOXLZ1136:         USE COMMONS, ONLY: NATOMS, NSAVE, PERIODIC, BOXLX, BOXLY, BOXLZ
1480: 1137: 
1481:         ! QMIN: energies of the saved minima1138:         ! QMIN: energies of the saved minima
1482:         ! QMINP: coordinates of the saved minima1139:         ! QMINP: coordinates of the saved minima
1483:         ! FF: step at which the saved minima was first found1140:         ! FF: step at which the saved minima was first found
1484:         USE QMODULE, ONLY: QMIN, QMINP, FF1141:         USE QMODULE, ONLY: QMIN, QMINP, FF
1485: 1142: 
1486:         ! BOX_DERIV: type for storing box derivative information 
1487:         ! CALC_BOX_DERIV: calculates box derivative information 
1488:         USE TRICLINIC_MOD, ONLY: BOX_DERIV, CALC_BOX_DERIV 
1489:  
1490:         IMPLICIT NONE1143:         IMPLICIT NONE
1491: 1144: 
1492:         ! GETUNIT: function for fiding a free file unit1145:         ! GETUNIT: function for fiding a free file unit
1493:         ! FUNIT: output file unit1146:         ! FUNIT: output file unit
1494:         ! NMOL: actual number of particles1147:         ! NMOL: actual number of particles
1495:         ! I, J: loop iterators1148:         ! I, J: loop iterators
1496:         INTEGER :: GETUNIT, FUNIT, NMOL, I, J1149:         INTEGER :: GETUNIT, FUNIT, NMOL, I, J
1497: 1150: 
1498:         ! BOX_SIZE: convenience array for the box sizes1151:         ! BOX_SIZE: convenience array for the box sizes
1499:         ! BOX_ANGLES: convenience array for the box angles1152:         ! BOX_ANGLES: convenience array for the box angles


r33077/potential.f90 2017-07-26 19:30:29.188736660 +0100 r33076/potential.f90 2017-07-26 19:30:31.360765503 +0100
 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 49:    USE OPP_MOD, ONLY: OPP
  50:    USE LJ_VAR_RAD_MOD, ONLY: LJ_VAR_RAD
 50:    USE EWALD 51:    USE EWALD
 51:    USE watermethane_mod 52:    USE watermethane_mod
 52:    IMPLICIT NONE 53:    IMPLICIT NONE
 53: ! Arguments 54: ! Arguments
 54: ! TODO: Work out intents 55: ! TODO: Work out intents
 55: ! TODO: Fix array dimensions? 56: ! TODO: Fix array dimensions?
 56:    REAL(REAL64)               :: X(*) 57:    REAL(REAL64)               :: X(*)
 57:    REAL(REAL64)               :: GRAD(*) 58:    REAL(REAL64)               :: GRAD(*)
 58:    REAL(REAL64)               :: EREAL 59:    REAL(REAL64)               :: EREAL
 59:    LOGICAL, INTENT(IN)        :: GRADT 60:    LOGICAL, INTENT(IN)        :: GRADT
1055: ! ds656> Multicomponent LJ with no cutoff (akin to BLJ_CLUST).1056: ! ds656> Multicomponent LJ with no cutoff (akin to BLJ_CLUST).
1056: !        Should give same results as GLJ, but different code.1057: !        Should give same results as GLJ, but different code.
1057:       CALL RAD(X, GRAD, EREAL, GRADT)1058:       CALL RAD(X, GRAD, EREAL, GRADT)
1058:       CALL MLJ(X, GRAD, EREAL, GRADT, SECT, .FALSE.)1059:       CALL MLJ(X, GRAD, EREAL, GRADT, SECT, .FALSE.)
1059: 1060: 
1060:    ELSE IF (GLJT) THEN1061:    ELSE IF (GLJT) THEN
1061: ! generalised LJ with no cutoff1062: ! generalised LJ with no cutoff
1062:       CALL RAD(X, GRAD, EREAL, GRADT)1063:       CALL RAD(X, GRAD, EREAL, GRADT)
1063:       CALL GLJ(X, GRAD, EREAL, GRADT)1064:       CALL GLJ(X, GRAD, EREAL, GRADT)
1064: 1065: 
 1066:    ELSE IF (LJ_VAR_RADT) THEN
 1067: ! jwrm2> Lennard-Jones with each particle having a different radius. Works with
 1068: !        RIGIDINIT
 1069:       IF (RIGIDINIT .AND. .NOT. ATOMRIGIDCOORDT) THEN
 1070:          XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)
 1071:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)
 1072:       END IF
 1073:       CALL LJ_VAR_RAD(X, GRAD, EREAL, GRADT)
 1074:       IF (RIGIDINIT .AND. .NOT. ATOMRIGIDCOORDT) THEN
 1075:          X(DEGFREEDOMS+1:3*NATOMS)=0.0D0
 1076:          X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS)
 1077:          IF (GRADT) THEN
 1078:             GRADATOMS(1:3*NATOMS) = GRAD(1:3*NATOMS) ! Needed for AACONVERGENCE
 1079:             CALL TRANSFORMGRAD(GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD)
 1080:             GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0
 1081:             GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)
 1082:          END IF
 1083:      END IF
 1084: 
1065:    ELSE IF (BINARY) THEN1085:    ELSE IF (BINARY) THEN
1066:       IF (SHIFTCUT) THEN1086:       IF (SHIFTCUT) THEN
1067:          if (rigidinit.and.(.not.atomrigidcoordt)) then1087:          if (rigidinit.and.(.not.atomrigidcoordt)) then
1068:             xrigidcoords(1:degfreedoms) = x(1:degfreedoms)1088:             xrigidcoords(1:degfreedoms) = x(1:degfreedoms)
1069:             call transformrigidtoc(1, nrigidbody, x, xrigidcoords)1089:             call transformrigidtoc(1, nrigidbody, x, xrigidcoords)
1070:          endif1090:          endif
1071:          call ljpshift(x, grad, ereal, gradt, sect)1091:          call ljpshift(x, grad, ereal, gradt, sect)
1072:          if (rigidinit.and.(.not.atomrigidcoordt)) then1092:          if (rigidinit.and.(.not.atomrigidcoordt)) then
1073:             x(degfreedoms+1:3*natoms)=0.0d01093:             x(degfreedoms+1:3*natoms)=0.0d0
1074:             x(1:degfreedoms)=xrigidcoords(1:degfreedoms)1094:             x(1:degfreedoms)=xrigidcoords(1:degfreedoms)
1091:       CALL LB2(X, GRAD, EREAL, GRADT)1111:       CALL LB2(X, GRAD, EREAL, GRADT)
1092: 1112: 
1093:    ELSE IF (LJCOULT) THEN1113:    ELSE IF (LJCOULT) THEN
1094:       CALL RADR(X, GRAD, EREAL, GRADT)1114:       CALL RADR(X, GRAD, EREAL, GRADT)
1095:       CALL LJCOUL(X, GRAD, EREAL, GRADT)1115:       CALL LJCOUL(X, GRAD, EREAL, GRADT)
1096: 1116: 
1097:   ELSE IF (LJ_GAUSST) THEN1117:   ELSE IF (LJ_GAUSST) THEN
1098:       CALL LJ_GAUSS(X, GRAD, EREAL, GRADT)1118:       CALL LJ_GAUSS(X, GRAD, EREAL, GRADT)
1099: 1119: 
1100:   ELSE IF (OPPT) THEN1120:   ELSE IF (OPPT) THEN
1101:       CALL OPP(X, GRAD, EREAL, GRADT, SECT)1121:       CALL OPP(X, GRAD, EREAL, GRADT)
1102: 1122: 
1103:    ELSE IF (STOCKT) THEN1123:    ELSE IF (STOCKT) THEN
1104:       CALL RADR(X, GRAD, EREAL, GRADT)1124:       CALL RADR(X, GRAD, EREAL, GRADT)
1105:       CALL STOCK(X, GRAD, EREAL, GRADT)1125:       CALL STOCK(X, GRAD, EREAL, GRADT)
1106: 1126: 
1107: ! RBAA potentials1127: ! RBAA potentials
1108: 1128: 
1109:    ELSE IF (CAPBINT) THEN1129:    ELSE IF (CAPBINT) THEN
1110:       CALL CAPBIN (X, GRAD, EREAL, GRADT)1130:       CALL CAPBIN (X, GRAD, EREAL, GRADT)
1111: 1131: 


r33077/triclinic.f90 2017-07-26 19:30:29.420739742 +0100 r33076/triclinic.f90 2017-07-26 19:30:31.588768532 +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/triclinic.f90' in revision 33076
  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 to calculatye transformation matrices and  
 21: ! derivatives for a triclinic unit cell, as well as a type to store them. 
 22: ! Questions to jwrm2 
 23: ! 
 24: MODULE TRICLINIC_MOD 
 25:  
 26: PUBLIC :: BOX_DERIV, CALC_BOX_DERIV 
 27:  
 28: ! BOX_DERIV: stores information for calculating lattice derivatives 
 29: TYPE :: BOX_DERIV 
 30:     ! ANGLES: unit cell angles 
 31:     ! SIZES: unit cell lengths 
 32:     DOUBLE PRECISION :: ANGLES(3), SIZES(3) 
 33:  
 34:     ! ORTHOG: the orthogonal tansformation matrix, for changing from 
 35:     !         crystallographic space to orthogonal space first index is row, 
 36:     !         second is column 
 37:     ! ORTHOGSQ: MATMUL(ORTHOG, ORTHOG), needed for second derivatives 
 38:     DOUBLE PRECISION :: ORTHOG(3,3), ORTHOGSQ(3, 3) 
 39:     ! ORTHOG: the inverse orthogonal tansformation matrix, for changing from 
 40:     !             orthogonal space to crystallographic space 
 41:     DOUBLE PRECISION :: ORTHOG_INV(3,3) 
 42:  
 43:     ! DERIV: derivatives of the orthogonal matrix wrt lattice parameters, 
 44:     !        first index indicates the row, second index the column and 
 45:     !        third index the parameter, with 1-3 being the angles and 4-6 being 
 46:     !        the lengths 
 47:     ! DERIV2: derivatives of DERIV wrt the six unit cell parameters 
 48:     DOUBLE PRECISION :: DERIV(3, 3, 6), DERIV2(3, 3, 6, 6) 
 49:  
 50:     ! DERIV_INV: derivatives of the inverse orthogonal matrix wrt lattice  
 51:     !        parameters, first index indicates the row, second index the column 
 52:     !        and third index the parameter, with 1-3 being the angles and 4-6  
 53:     !        being the lengths 
 54:     DOUBLE PRECISION :: DERIV_INV(3, 3, 6) 
 55:  
 56:     ! RECIP: lengths of the reciprocal lattice vectors 
 57:     DOUBLE PRECISION :: RECIP(3) 
 58:  
 59:     ! CELL_RANGE: number of cells in each direction that must be evaluated 
 60:     INTEGER :: CELL_RANGE(3) 
 61:  
 62:     ! V: dimensionless volume of the unit cell, ie volume if all lengths were 1 
 63:     DOUBLE PRECISION :: V 
 64:  
 65:     ! V_DERIV: derivatives of V wrt the lattice angles 
 66:     ! V_DERIV2: sedond derivatives of V wrt the lattice angles 
 67:     DOUBLE PRECISION :: V_DERIV(3), V_DERIV2(3, 3) 
 68: END TYPE BOX_DERIV 
 69:  
 70: CONTAINS 
 71:  
 72: !------------------------------------------------------------------------------- 
 73: ! 
 74: ! Calculates factors for the lattice derivatives that only depend on the cell 
 75: ! parameters. 
 76: ! ANGLES: alpha, beta and gamma lattice parameters 
 77: ! SIZES: a, b and c lattice parameters 
 78: ! CUTOFF: (optional) provide for working out the number of periodic repeats in 
 79: !         each direction within the cutoff sphere 
 80: ! 
 81: !------------------------------------------------------------------------------- 
 82:     PURE TYPE(BOX_DERIV) FUNCTION CALC_BOX_DERIV(ANGLES, SIZES, CUTOFF) 
 83:  
 84:         IMPLICIT NONE 
 85:  
 86:         DOUBLE PRECISION, INTENT(IN) :: ANGLES(3), SIZES(3) 
 87:         DOUBLE PRECISION, OPTIONAL, INTENT(IN) :: CUTOFF 
 88:  
 89:         ! V: factor, related to (but not quite equal to) the unit cell volume 
 90:         ! DVDA: derivatives of V wrt the three lattice angles 
 91:         ! C, S: cos and sin of the lattice angles 
 92:         ! AL, BL, CL: unit cell lengths 
 93:         ! ORTHOG, INV_ORTHOG, DERIV, DERIV2, DERIV_INV: parts of the output 
 94:         DOUBLE PRECISION :: V, DVDA(3), C(3), S(3), AL, BL, CL 
 95:         DOUBLE PRECISION :: ORTHOG(3, 3), ORTHOG_INV(3, 3) 
 96:         DOUBLE PRECISION :: DERIV(3, 3, 6), DERIV2(3,3,6,6), DERIV_INV(3, 3, 6) 
 97:  
 98:         ! Initialise output variables 
 99:         ORTHOG(:,:) = 0.D0 
100:         DERIV(:,:,:) = 0.D0 
101:         DERIV2(:,:,:,:) = 0.D0 
102:         ORTHOG_INV(:,:) = 0.D0 
103:         DERIV_INV(:,:,:) = 0.D0 
104:  
105:         ! Calculate factors 
106:         AL = SIZES(1) 
107:         BL = SIZES(2) 
108:         CL = SIZES(3) 
109:         C(:) = COS(ANGLES(:)) 
110:         S(:) = SIN(ANGLES(:)) 
111:         V = SQRT(1 - C(1)**2 - C(2)**2 - C(3)**2 + 2.D0 * C(1) * C(2) * C(3)) 
112:         DVDA(1) = S(1) * (C(1) - C(2) * C(3)) / V 
113:         DVDA(2) = S(2) * (C(2) - C(3) * C(1)) / V 
114:         DVDA(3) = S(3) * (C(3) - C(1) * C(2)) / V 
115:  
116:         ! Calculate the orthogonal transformation matrix 
117:         ORTHOG(1, 1) = AL 
118:         ORTHOG(1, 2) = BL * C(3) 
119:         ORTHOG(1, 3) = CL * C(2) 
120:         ORTHOG(2, 2) = BL * S(3) 
121:         ORTHOG(2, 3) = CL * (C(1) - C(2) * C(3)) / S(3) 
122:         ORTHOG(3, 3) = CL * V / S(3) 
123:  
124:         ! Calculate the derivatives of the othogonal matrix 
125:         ! Derivatives of the top row wrt angles 
126:         DERIV(1, 3, 2) = - CL * S(2) 
127:         DERIV(1, 2, 3) = - BL * S(3) 
128:         ! Derivatives of the top row wrt to lengths 
129:         DERIV(1, 1, 4) = 1.D0 
130:         DERIV(1, 2, 5) = C(3) 
131:         DERIV(1, 3, 6) = C(2) 
132:         ! Derivatives of the middle row wrt angles 
133:         DERIV(2, 3, 1) = - CL * S(1) / S(3) 
134:         DERIV(2, 3, 2) = CL * S(2) * C(3) / S(3) 
135:         DERIV(2, 2, 3) = BL * C(3) 
136:         DERIV(2, 3, 3) = CL * (C(2) - C(1) * C(3)) / S(3)**2 
137:         ! Derivatives of the middle row wrt to lengths 
138:         DERIV(2, 2, 5) = S(3) 
139:         DERIV(2, 3, 6) = (C(1) - C(2) * C(3)) / S(3) 
140:         ! Derivatives of the bottom row wrt angles 
141:         DERIV(3, 3, 1) = CL * DVDA(1) / S(3) 
142:         DERIV(3, 3, 2) = CL * DVDA(2) / S(3) 
143:         DERIV(3, 3, 3) = CL * (DVDA(3) - C(3) * V / S(3)) / S(3) 
144:         ! Derivatives of the bottom row wrt lengths 
145:         DERIV(3, 3, 6) = V / S(3) 
146:  
147:         ! Calculate the second derivatives of the orthogonal matrix. FML. 
148:         ! Derivatives of element 12 
149:         DERIV2(1, 2, 3, 3) = - BL * C(3) 
150:         DERIV2(1, 2, 3, 5) = - S(3) 
151:         DERIV2(1, 2, 5, 3) = DERIV2(1, 2, 3, 5) 
152:         ! Derivatives of element 13 
153:         DERIV2(1, 3, 2, 2) = - CL * C(2) 
154:         DERIV2(1, 3, 2, 6) = - S(2) 
155:         DERIV2(1, 3, 6, 2) = DERIV2(1, 3, 2, 6) 
156:         ! Derivatives of element 22 
157:         DERIV2(2, 2, 3, 3) = - BL * S(3) 
158:         DERIV2(2, 2, 3, 5) = C(3) 
159:         DERIV2(2, 2, 5, 3) = DERIV2(2, 2, 3, 5) 
160:         ! Derivatives of element 23 
161:         DERIV2(2, 3, 1, 1) = - CL * C(1) / S(3) 
162:         DERIV2(2, 3, 1, 3) = CL * S(1) * C(3) / S(3)**2 
163:         DERIV2(2, 3, 3, 1) = DERIV2(2, 3, 1, 3) 
164:         DERIV2(2, 3, 1, 6) = - S(1) / S(3) 
165:         DERIV2(2, 3, 6, 1) = DERIV2(2, 3, 1, 6) 
166:         DERIV2(2, 3, 2, 2) = CL * C(2) * C(3) / S(3) 
167:         DERIV2(2, 3, 2, 3) = - CL * S(2) / S(3)**2 
168:         DERIV2(2, 3, 3, 2) = DERIV2(2, 3, 2, 3) 
169:         DERIV2(2, 3, 2, 6) = S(2) * C(3) / S(3) 
170:         DERIV2(2, 3, 6, 2) = DERIV2(2, 3, 2, 6) 
171:         DERIV2(2, 3, 3, 3) = CL*(2.D0 * C(3) * (C(1) * C(3) - C(2)) / S(3)**3 +& 
172:             C(1) / S(3)) 
173:         DERIV2(2, 3, 3, 6) = (C(2) - C(1) * C(3)) / S(3)**2 
174:         DERIV2(2, 3, 6, 3) = DERIV2(2, 3, 3, 6) 
175:         ! Derivatives of element 33 
176:         DERIV2(3, 3, 1, 1) = CL * ((C(1)**2 - S(1)**2 - C(1) * C(2) * C(3)) - & 
177:             DVDA(1)**2) / (V * S(3)) 
178:         DERIV2(3, 3, 1, 2) = CL * (S(1) * S(2) * C(3) - DVDA(1) * DVDA(2)) / & 
179:             (V * S(3)) 
180:         DERIV2(3, 3, 2, 1) = DERIV2(3, 3, 1, 2) 
181:         DERIV2(3, 3, 1, 3) = CL * (S(1) * C(2) / V - & 
182:             DVDA(1) * DVDA(3) / (S(3) * V) - DVDA(1) * C(3) / S(3)**2) 
183:         DERIV2(3, 3, 3, 1) = DERIV2(3, 3, 1, 3) 
184:         DERIV2(3, 3, 1, 6) = DVDA(1) / S(3) 
185:         DERIV2(3, 3, 6, 1) = DERIV2(3, 3, 1, 6) 
186:         DERIV2(3, 3, 2, 2) = CL * ((C(2)**2 - S(2)**2 - C(1) * C(2) * C(3)) - & 
187:             DVDA(2)**2) / (V * S(3)) 
188:         DERIV2(3, 3, 2, 3) = CL * (S(2) * C(1) / V - & 
189:             DVDA(2) * DVDA(3) / (S(3) * V) - DVDA(2) * C(3) / S(3)**2) 
190:         DERIV2(3, 3, 3, 2) = DERIV2(3, 3, 2, 3) 
191:         DERIV2(3, 3, 2, 6) = DVDA(2) / S(3) 
192:         DERIV2(3, 3, 6, 2) = DERIV2(3, 3, 2, 6) 
193:         DERIV2(3, 3, 3, 3) = CL * ((C(3)**2 - S(3)**2 - C(1) * C(2) * C(3))/V -& 
194:             DVDA(3)**2 / V - C(3) * DVDA(3) / S(3) + V / S(3)**2)/ S(3) - & 
195:             CL * C(3) * (DVDA(3) - V * C(3) / S(3)) / S(3)**2 
196:         DERIV2(3, 3, 3, 6) = (DVDA(3) - V * C(3) / S(3)) / S(3) 
197:         DERIV2(3, 3, 6, 3) = DERIV2(3, 3, 3, 6) 
198:  
199:         ! The inverse of the orthogonal transformation matrix 
200:         ORTHOG_INV(1, 1) = 1.D0 / AL 
201:         ORTHOG_INV(1, 2) = - C(3) / ( AL * S(3)) 
202:         ORTHOG_INV(1, 3) = (C(1) * C(3) - C(2)) / (AL * S(3) * V) 
203:         ORTHOG_INV(2, 2) = 1.D0 / (BL * S(3)) 
204:         ORTHOG_INV(2, 3) = (C(2) * C(3) - C(1)) / (BL * S(3) * V) 
205:         ORTHOG_INV(3, 3) = S(3) / (CL * V) 
206:  
207:         ! Calculate the derivatives of the inverse othogonal matrix 
208:         ! Derivatives of the top row wrt angles 
209:         DERIV_INV(1, 3, 1) = - S(1) * C(3) / (AL * S(3) * V) + & 
210:             (C(1) * C(3) - C(2)) * (C(2) * C(3) - C(1))*S(1)/(AL * S(3) * V**3) 
211:         DERIV_INV(1, 3, 2) = S(2) / (AL * S(3) * V) + & 
212:             (C(1) * C(3) - C(2))**2 * S(2) / (AL * S(3) * V**3) 
213:         DERIV_INV(1, 2, 3) = 1.D0 / (AL * S(3)**2) 
214:         DERIV_INV(1, 3, 3) = - C(1) / (AL * V) + & 
215:             (C(1) * C(3) - C(2)) * ((C(1) * C(2) - C(3)) / & 
216:             (AL * V**3) - C(3) / (AL * S(3)**2 * V)) 
217:         ! Derivatives of the top row wrt to lengths 
218:         DERIV_INV(1, 1, 4) = - 1.D0 / AL**2 
219:         DERIV_INV(1, 2, 4) = C(3) / (AL**2 * S(3)) 
220:         DERIV_INV(1, 3, 4) = (C(2) - C(1) * C(3)) / (AL**2 * S(3) * V) 
221:         ! Derivatives of the middle row wrt angles 
222:         DERIV_INV(2, 3, 1) = S(1) / (BL * S(3) * V) + & 
223:             (C(2) * C(3) - C(1))**2 * S(1) / (BL * S(3) * V**3) 
224:         DERIV_INV(2, 3, 2) = - S(2) * C(3) / (BL * S(3) * V) + & 
225:             (C(2) * C(3) - C(1)) * (C(1) * C(3) - C(2))*S(2)/(BL * S(3) * V**3) 
226:         DERIV_INV(2, 2, 3) = - C(3) / (BL * S(3)**2) 
227:         DERIV_INV(2, 3, 3) = - C(2) / (BL * V) + & 
228:             (C(2) * C(3) - C(1)) * ((C(1) * C(2) - C(3)) / & 
229:             (BL * V**3) - C(3) / (BL * S(3)**2 * V)) 
230:         ! Derivatives of the middle row wrt to lengths 
231:         DERIV_INV(2, 2, 5) = - 1.D0 / (BL**2 * S(3)) 
232:         DERIV_INV(2, 3, 5) = (C(1) - C(2) * C(3)) / (BL**2 * S(3) * V) 
233:         ! Derivatives of the bottom row wrt angles 
234:         DERIV_INV(3, 3, 1) = S(1) * S(3) * (C(2) * C(3) - C(1)) / (CL * V**3) 
235:         DERIV_INV(3, 3, 2) = S(2) * S(3) * (C(1) * C(3) - C(2)) / (CL * V**3) 
236:         DERIV_INV(3, 3, 3) = C(3) / (CL * V) + & 
237:             S(3)**2 * (C(1) * C(2) - C(3)) / (CL * V**3) 
238:         ! Derivatives of the bottom row wrt lengths 
239:         DERIV_INV(3, 3, 6) = - S(3) / (CL**2 * V) 
240:  
241:         ! Calculate the derivative of volume wrt the angles 
242:         CALC_BOX_DERIV%V_DERIV(1) = S(1) * (C(1) - C(2) * C(3)) / V 
243:         CALC_BOX_DERIV%V_DERIV(2) = S(2) * (C(2) - C(3) * C(1)) / V 
244:         CALC_BOX_DERIV%V_DERIV(3) = S(3) * (C(3) - C(1) * C(2)) / V 
245:  
246:         ! Calculate second derivatives of volume wrt the angles 
247:         CALC_BOX_DERIV%V_DERIV2(1,1) = (S(1) * (C(1) - C(2) * C(3))**2 - & 
248:             C(1)**2 + S(1)**2 - C(1) * C(2) * C(3)) / V**3 
249:         CALC_BOX_DERIV%V_DERIV2(2,2) = (S(2) * (C(2) - C(3) * C(1))**2 - & 
250:             C(2)**2 + S(2)**2 - C(2) * C(3) * C(1)) / V**3 
251:         CALC_BOX_DERIV%V_DERIV2(3,3) = (S(3) * (C(3) - C(1) * C(2))**2 - & 
252:             C(3)**2 + S(3)**2 - C(3) * C(1) * C(2)) / V**3 
253:         CALC_BOX_DERIV%V_DERIV2(1,2) = (C(1) * C(2) - & 
254:             C(3)) * S(1) * S(2) * S(3)**2 / V**3 
255:         CALC_BOX_DERIV%V_DERIV2(2,1) = CALC_BOX_DERIV%V_DERIV2(1,2) 
256:         CALC_BOX_DERIV%V_DERIV2(2,3) = (C(2) * C(3) - & 
257:             C(1)) * S(2) * S(3) * S(1)**2 / V**3 
258:         CALC_BOX_DERIV%V_DERIV2(3,2) = CALC_BOX_DERIV%V_DERIV2(2,3) 
259:         CALC_BOX_DERIV%V_DERIV2(3,1) = (C(3) * C(1) - & 
260:             C(2)) * S(3) * S(1) * S(2)**2 / V**3 
261:         CALC_BOX_DERIV%V_DERIV2(1,3) = CALC_BOX_DERIV%V_DERIV2(3,1) 
262:  
263:         ! Calculate the lengths of the reciprocal lattice vectors 
264:         CALC_BOX_DERIV%RECIP(:) = S(:) / (SIZES(:) * V) 
265:  
266:         ! Calculate the number of cells to check in each direction 
267:         IF (PRESENT(CUTOFF)) THEN 
268:             CALC_BOX_DERIV%CELL_RANGE(:) = CEILING(CUTOFF * & 
269:                 CALC_BOX_DERIV%RECIP(:)) 
270:         ELSE 
271:             CALC_BOX_DERIV%CELL_RANGE(:) = 0.D0 
272:         END IF 
273:  
274:         ! Set the output 
275:         CALC_BOX_DERIV%ANGLES(:) = ANGLES(:) 
276:         CALC_BOX_DERIV%SIZES(:) = SIZES(:) 
277:         CALC_BOX_DERIV%ORTHOG(:,:) = ORTHOG(:,:) 
278:         CALC_BOX_DERIV%ORTHOGSQ(:,:) = & 
279:             MATMUL(TRANSPOSE(ORTHOG(:,:)), ORTHOG(:,:)) 
280:         CALC_BOX_DERIV%DERIV(:,:,:) = DERIV(:,:,:) 
281:         CALC_BOX_DERIV%DERIV2(:,:,:,:) = DERIV2(:,:,:,:) 
282:         CALC_BOX_DERIV%ORTHOG_INV(:,:) = ORTHOG_INV(:,:) 
283:         CALC_BOX_DERIV%DERIV_INV(:,:,:) = DERIV_INV(:,:,:) 
284:         CALC_BOX_DERIV%V = V 
285:  
286:     END FUNCTION CALC_BOX_DERIV 
287:  
288: END MODULE TRICLINIC_MOD 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0