hdiff output

r32371/keywords.f 2017-04-25 16:30:09.842726889 +0100 r32370/keywords.f 2017-04-25 16:30:10.526735735 +0100
 34:      &                      rotcentre, rotpselect, rotoccuw, rotcutoff, setchiralgeneric, PRMTOP, IGB, RGBMAX, CUT, 34:      &                      rotcentre, rotpselect, rotoccuw, rotcutoff, setchiralgeneric, PRMTOP, IGB, RGBMAX, CUT,
 35:      &                      SALTCON, macroiont, nmacroions, macroiondist 35:      &                      SALTCON, macroiont, nmacroions, macroiondist
 36:       USE modamber 36:       USE modamber
 37:       USE PORFUNCS 37:       USE PORFUNCS
 38:       USE MYGA_PARAMS 38:       USE MYGA_PARAMS
 39:       USE BGUPMOD 39:       USE BGUPMOD
 40:       USE GLJYMOD 40:       USE GLJYMOD
 41:       USE CHIRO_MODULE, ONLY: CHIRO_SIGMA, CHIRO_MU, CHIRO_GAMMA, CHIRO_L 41:       USE CHIRO_MODULE, ONLY: CHIRO_SIGMA, CHIRO_MU, CHIRO_GAMMA, CHIRO_L
 42:       USE CONVEX_POLYHEDRA_MODULE, ONLY: INITIALISE_POLYHEDRA, K_COMPRESS, K_OVERLAP 42:       USE CONVEX_POLYHEDRA_MODULE, ONLY: INITIALISE_POLYHEDRA, K_COMPRESS, K_OVERLAP
 43:       USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS_MODE, LJ_GAUSS_RCUT, LJ_GAUSS_EPS, 43:       USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS_MODE, LJ_GAUSS_RCUT, LJ_GAUSS_EPS,
 44:      &                        LJ_GAUSS_R0, LJ_GAUSS_SIGMASQ, LJ_GAUSS_PARAMS, 44:      &                        LJ_GAUSS_R0, LJ_GAUSS_INITIALISE
 45:      &                        LJ_GAUSS_INITIALISE 
 46:       USE MBPOLMOD, ONLY: MBPOLINIT 45:       USE MBPOLMOD, ONLY: MBPOLINIT
 47:       USE SWMOD, ONLY: SWINIT, MWINIT 46:       USE SWMOD, ONLY: SWINIT, MWINIT
 48:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_GET_COORDS 47:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_GET_COORDS
 49: !     &                                 AMBER12_ATOMSS, AMBER12_SETUP, 48: !     &                                 AMBER12_ATOMSS, AMBER12_SETUP,
 50: !     &                                 AMBER12_RESIDUES, 49: !     &                                 AMBER12_RESIDUES,
 51: !     &                                 POPULATE_ATOM_DATA 50: !     &                                 POPULATE_ATOM_DATA
 52:       USE CHIRALITY, ONLY : CIS_TRANS_TOL 51:       USE CHIRALITY, ONLY : CIS_TRANS_TOL
 53:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR 52:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR
 54:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS, 53:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS,
 55:      &     PARSE_MLJ_PARAMS 54:      &     PARSE_MLJ_PARAMS
6227:          ALLOCATE(SITE(NRBSITES,3))6226:          ALLOCATE(SITE(NRBSITES,3))
6228: !        Maybe the above two lines are not necessary!6227: !        Maybe the above two lines are not necessary!
6229: 6228: 
6230: ! jwrm2>6229: ! jwrm2>
6231:       ELSE IF (WORD.EQ.'LJGAUSS') THEN6230:       ELSE IF (WORD.EQ.'LJGAUSS') THEN
6232:           LJ_GAUSST = .TRUE.6231:           LJ_GAUSST = .TRUE.
6233:           CALL READI(LJ_GAUSS_MODE)6232:           CALL READI(LJ_GAUSS_MODE)
6234:           CALL READF(LJ_GAUSS_RCUT)6233:           CALL READF(LJ_GAUSS_RCUT)
6235:           CALL READF(LJ_GAUSS_EPS)6234:           CALL READF(LJ_GAUSS_EPS)
6236:           CALL READF(LJ_GAUSS_R0)6235:           CALL READF(LJ_GAUSS_R0)
6237:           CALL READF(LJ_GAUSS_SIGMASQ) 
6238:           IF (LJ_GAUSS_MODE .EQ. 3) CALL READI(LJ_GAUSS_PARAMS) 
6239:           CALL LJ_GAUSS_INITIALISE()6236:           CALL LJ_GAUSS_INITIALISE()
6240: ! jwrm2> end6237: ! jwrm2> end
6241: 6238: 
6242:       ELSE IF (WORD.EQ.'STOCK') THEN6239:       ELSE IF (WORD.EQ.'STOCK') THEN
6243:          STOCKT=.TRUE.6240:          STOCKT=.TRUE.
6244:          RIGID=.TRUE.6241:          RIGID=.TRUE.
6245:          NRBSITES=16242:          NRBSITES=1
6246:          CALL READF(STOCKMU)6243:          CALL READF(STOCKMU)
6247:          CALL READF(STOCKLAMBDA)6244:          CALL READF(STOCKLAMBDA)
6248:          ALLOCATE(SITE(NRBSITES,3))6245:          ALLOCATE(SITE(NRBSITES,3))


r32371/lj_gauss.f90 2017-04-25 16:30:10.066729785 +0100 r32370/lj_gauss.f90 2017-04-25 16:30:10.750738631 +0100
 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of
 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13: !   GNU General Public License for more details. 13: !   GNU General Public License for more details.
 14: ! 14: !
 15: !   You should have received a copy of the GNU General Public License 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 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 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: ! 18: !
 19: ! 19: !
 20: ! This module provides a routine for calculating the energy and gradients with 20: ! This module provides a routine for calculating the energy and gradients with
 21: ! a double well Lennard-Jones + Gaussian potential. The possibility of varying 21: ! a double well Lennard-Jones + Gaussian potential.
 22: ! and optimising the potential parameneter is included. 
 23: ! Also included is a framework for a varying triclinic unit cell. This could be 
 24: ! separated out and generalised to all isotropic pair potentials. 
 25: ! Questions to jwrm2. 22: ! Questions to jwrm2.
 26: ! 23: !
 27: MODULE LJ_GAUSS_MOD 24: MODULE LJ_GAUSS_MOD
 28:  25: 
 29: ! Public parameters 26: ! Public parameters
 30: PUBLIC :: LJ_GAUSS_MODE, LJ_GAUSS_RCUT, LJ_GAUSS_EPS, LJ_GAUSS_R0 27: PUBLIC :: LJ_GAUSS_MODE, LJ_GAUSS_RCUT, LJ_GAUSS_EPS, LJ_GAUSS_R0
 31: PUBLIC :: LJ_GAUSS_SIGMASQ, LJ_GAUSS_PARAMS 
 32: ! Public functions 28: ! Public functions
 33: PUBLIC :: LJ_GAUSS, LJ_GAUSS_INITIALISE, VIEW_LJ_GAUSS, LJ_GAUSS_TAKESTEP 29: PUBLIC :: LJ_GAUSS, LJ_GAUSS_INITIALISE, VIEW_LJ_GAUSS
 34: ! Private parameters 30: ! Private parameters
 35: PRIVATE :: RCUTINV6, RCUTINV12, EXP_RCUT, PREF_RCUT, PI, SCL_FCT 31: PRIVATE :: LJ_GAUSS_SIGMASQ, RCUTINV6, RCUTINV12, EXP_RCUT, PREF_RCUT, PI
 36: PRIVATE :: IMAGE_CUTOFF, MAX_LENGTH_STEP, MAX_ANGLE_STEP, MAX_EPS_STEP 32: PRIVATE :: SCL_FCT
 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 
 39: PRIVATE :: SIGMASQ_UPPER, SIGMASQ_REPULSION 
 40: ! Private functions 33: ! Private functions
 41: PRIVATE :: PERIODIC_RESET, LJ_GAUSS_PER, CALC_ENERGY, CALC_ENERGY_LJ 34: PRIVATE :: PERIODIC_RESET, LJ_GAUSS_PER, CALC_ENERGY, CALC_ENERGY_LJ
 42: PRIVATE :: CALC_ENERGY_GAUSS, CALC_GRAD, CALC_FCTS, CONSTRAIN_PARAMETERS 35: PRIVATE :: CALC_ENERGY_GAUSS, CALC_GRAD, CALC_FCTS, CONSTRAIN_PARAMETERS
 43: PRIVATE :: BOX_DERIV, CALC_BOX_DERIV, CHECK_ANGLES, REJECT, CONSTRAIN_VOLUME 36: PRIVATE :: BOX_DERIV, CALC_BOX_DERIV, CHECK_ANGLES, REJECT
 44:  37: 
 45: ! EPS: strength of the Gaussian potential. The LJ strength is 1. 38: ! EPS: strength of the Gaussian potential. The LJ strength is 1.
 46: ! R0: Equilibrium distance of the minimum in the Gaussian potential. The LJ 39: ! R0: Equilibrium distance of the minimum in the Gaussian potential. The LJ
 47: !     rmin is 1. 40: !     rmin is 1.
 48: ! SIGMASQ: variation (width) of the Gaussian potential 
 49: ! RCUTSQ: cutoff distance at which the energy, gradients and Hessian of 41: ! RCUTSQ: cutoff distance at which the energy, gradients and Hessian of
 50: !         both parts of the potential go smoothly to zero. 42: !         both parts of the potential go smoothly to zero.
 51: DOUBLE PRECISION :: LJ_GAUSS_EPS, LJ_GAUSS_R0, LJ_GAUSS_SIGMASQ, LJ_GAUSS_RCUT 43: DOUBLE PRECISION :: LJ_GAUSS_EPS, LJ_GAUSS_R0, LJ_GAUSS_RCUT
 52:  44: 
 53: ! LJ_GAUSS_MODE: 0 standard minimisation 45: ! LJ_GAUSS_MODE: 0 standard minimisation
 54: !                1 the final triplet of coordinates will be interpreted as the 46: !                1 the final triplet of coordinates will be interpreted as the
 55: !                  unit cell parameters of an orthogonal unit cell, and 47: !                  unit cell parameters of an orthogonal unit cell, and
 56: !                  optimised 48: !                  optimised
 57: !                2 the final triplet of coordinates are the unit cell length 49: !                2 the final triplet of coordinates are the unit cell length
 58: !                  parameters and the second last triplet are the unit cell 50: !                  parameters and the second last triplet are the unit cell
 59: !                  angles, giving a triclinic unit cell, which will be optimised 51: !                  angles, giving a triclinic unit cell, which will be optimised
 60: !                3 as for 2, but the third last triplet are the potential 52: !                3 as for 2, but the third last triplet are the potential
 61: !                  parameters LJ_GAUSS_EPS, LJ_GAUSS_R0, LJ_GAUSS_SIGMASQ and 53: !                  parameters LJ_GAUSS_EPS, LJ_GAUSS_R0, 0.D0 and the potential
 62: !                  the potential parameters will be optimised, depending on the 54: !                  parameters will be optimised
 63: !                  value of LJ_GAUSS_PARAMS 55: INTEGER :: LJ_GAUSS_MODE
 64: ! LJ_GAUSS_PARAMS: Specifies which of the three potential parameters will be 56: 
 65: !                  optimised. An integer between 1 and 7. 57: ! LJ_GAUSS_SIGMA: variation (width) of the Gaussian potential
 66: !                  1 bit set: LJ_GAUSS_EPS will be optimised 58: DOUBLE PRECISION, PARAMETER :: LJ_GAUSS_SIGMASQ = 0.02D0
 67: !                  2 bit set: LJ_GAUSS_R0 will be optimised 
 68: !                  3 bit set: LJ_GAUSS_SIGMASQ will be optimised 
 69: INTEGER :: LJ_GAUSS_MODE, LJ_GAUSS_PARAMS 
 70:  59: 
 71: ! RCUTINV6, RCUTINV12: factors for LJ which don't pened on particle distance 60: ! RCUTINV6, RCUTINV12: factors for LJ which don't pened on particle distance
 72: DOUBLE PRECISION :: RCUTINV6, RCUTINV12 61: DOUBLE PRECISION :: RCUTINV6, RCUTINV12
 73:  62: 
 74: ! EXP_RCUT, PREF_RCUT: factors for Gauss which don't depend on the particle 63: ! EXP_RCUT, PREF_RCUT: factors for Gauss which don't depend on the particle
 75: !                      distance 64: !                      distance
 76: ! SCL_FCT: overall scale factor for the potential, to allow comparison between 65: ! SCL_FCT: overall scale factor for the potential, to allow comparison between
 77: !          different parameters 66: !          different parameters
 78: DOUBLE PRECISION :: EXP_RCUT, PREF_RCUT, SCL_FCT 67: DOUBLE PRECISION :: EXP_RCUT, PREF_RCUT, SCL_FCT
 79:  68: 
 80: ! PI: the usual meaning 69: ! PI: the usual meaning
 81: DOUBLE PRECISION, PARAMETER :: PI = 4*ATAN(1.D0) 70: DOUBLE PRECISION, PARAMETER :: PI = 4*ATAN(1.D0)
 82:  71: 
 83: ! IMAGE_CUTOFF: if more periodic images than this in any direction fall 
 84: !               within the cutoff, the evaluation will not be carried 
 85: !               out and the step will be rejected. 
 86: DOUBLE PRECISION, PARAMETER :: IMAGE_CUTOFF = 30 
 87:  
 88: ! Parameters specifying the largest steps in LJ_GAUSS_TAKESTEP 
 89: ! MAX_LENGTH_STEP: largest allowed step in the length of a lattice 
 90: !                  vector 
 91: ! MAX_ANGLE_STEP: largest allowed step in the angle of a lattice vector 
 92: ! MAX_EPS_STEP: largest allowed step in LJ_GAUSS_EPS 
 93: ! MAX_R0_STEP: largest allowed step in LJ_GAUSS_R0 
 94: ! MAX_SIGMASQ_STEP: largest allowed step in LJ_GAUSS_SIGMASQ 
 95: DOUBLE PRECISION, PARAMETER :: MAX_LENGTH_STEP = 0.3D0 
 96: DOUBLE PRECISION, PARAMETER :: MAX_ANGLE_STEP = 0.1D0 
 97: DOUBLE PRECISION, PARAMETER :: MAX_EPS_STEP = 0.2D0 
 98: DOUBLE PRECISION, PARAMETER :: MAX_R0_STEP = 0.1D0 
 99: DOUBLE PRECISION, PARAMETER :: MAX_SIGMASQ_STEP = 0.001D0 
100:  
101: ! Parameters used for constraining the unit cell volume and the potential 
102: ! parameters 
103: ! V_EPS: scaling factor for the unit cell volume contraint energy 
104: ! V_SIGMA: WCA cutoff distance 
105: ! x_LOWER, x_UPPER: bounds on potential parameters 
106: ! x_REPULSION: harmonic force constant for each potential parameter 
107: DOUBLE PRECISION, PARAMETER :: V_EPS = 1.D-3 
108: DOUBLE PRECISION, PARAMETER :: V_SIGMA = 3.D-1 
109: DOUBLE PRECISION, PARAMETER :: EPS_LOWER = 0.D0, EPS_UPPER = 5.D0 
110: DOUBLE PRECISION, PARAMETER :: EPS_REPULSION = 1.D4 
111: DOUBLE PRECISION, PARAMETER :: R0_LOWER = 1.D0, R0_UPPER = 2.1D0 
112: DOUBLE PRECISION, PARAMETER :: R0_REPULSION = 1.D4 
113: DOUBLE PRECISION, PARAMETER :: SIGMASQ_LOWER = 0.01D0, SIGMASQ_UPPER = 0.03D0 
114: DOUBLE PRECISION, PARAMETER :: SIGMASQ_REPULSION = 1.D6 
115:  
116: ! BOX_DERIV: stores information for calculating lattice derivatives 72: ! BOX_DERIV: stores information for calculating lattice derivatives
117: TYPE :: BOX_DERIV 73: TYPE :: BOX_DERIV
118:     ! ORTHOG: the orthogonal tansformation matrix, for changing from 74:     ! ORTHOG: the orthogonal tansformation matrix, for changing from
119:     !         crystallographic space to orthogonal space first index is row, 75:     !         crystallographic space to orthogonal space first index is row,
120:     !         second is column 76:     !         second is column
121:     DOUBLE PRECISION :: ORTHOG(3,3) 77:     DOUBLE PRECISION :: ORTHOG(3,3)
122:  78: 
123:     ! DERIV: derivatives of the orthogonal matrix wrt lattice parameters, 79:     ! DERIV: derivatives of the orthogonal matrix wrt lattice parameters,
124:     !        first index indicates the row, second index the column and 80:     !        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 81:     !        third index the parameter, with 1-3 being the angles and 4-6 being
126:     !        the lengths 82:     !        the lengths
127:     DOUBLE PRECISION :: DERIV(3, 3, 6) 83:     DOUBLE PRECISION :: DERIV(3, 3, 6)
128:  84: 
129:     ! RECIP: lengths of the reciprocal lattice vectors 85:     ! RECIP: lengths of the reciprocal lattice vectors
130:     DOUBLE PRECISION :: RECIP(3) 86:     DOUBLE PRECISION :: RECIP(3)
131:  87: 
132:     ! CELL_RANGE: number of cells in each direction that must be evaluated 88:     ! CELL_RANGE: number of cells in each direction that must be evaluated
133:     INTEGER :: CELL_RANGE(3) 89:     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 90: END TYPE BOX_DERIV
141:  91: 
142: CONTAINS 92: CONTAINS
143:  93: 
144: !------------------------------------------------------------------------------- 94: !-------------------------------------------------------------------------------
145: ! 95: !
146: ! Main routine for calculating the energy and gradients. 96: ! Main routine for calculating the energy and gradients.
147: ! X: coordinates array 97: ! X: coordinates array
148: ! GRAD: gradients array 98: ! GRAD: gradients array
149: ! ENERGY: calculated energy 99: ! ENERGY: calculated energy
168:         INTEGER :: I, J, NMOL118:         INTEGER :: I, J, NMOL
169: 119: 
170:         ! RIJ: Interparticle vector120:         ! RIJ: Interparticle vector
171:         ! MODRIJ: distance between a pair of particles121:         ! MODRIJ: distance between a pair of particles
172:         ! DVDR: derivative of pair potential wrt MODRIJ122:         ! DVDR: derivative of pair potential wrt MODRIJ
173:         DOUBLE PRECISION :: RIJ(3), MODRIJ, DVDR123:         DOUBLE PRECISION :: RIJ(3), MODRIJ, DVDR
174: 124: 
175:         ! Factors for triclinic lattice125:         ! Factors for triclinic lattice
176:         TYPE(BOX_DERIV) :: BD126:         TYPE(BOX_DERIV) :: BD
177: 127: 
178: !        DO I = 1, NATOMS128:  !       DO I = 1, NATOMS
179: !            WRITE (MYUNIT, *) X(3*I-2:3*I)129:  !           WRITE (MYUNIT, *) X(3*I-2:3*I)
180: !        END DO130:  !       END DO
181: 131: 
182:         ! Initialise output variables132:         ! Initialise output variables
183:         ENERGY = 0.D0133:         ENERGY = 0.D0
184:         GRAD(:) = 0.D0134:         GRAD(:) = 0.D0
185: 135: 
186:         ! Calculate factors that do depend on parameters, but not on RIJ136:         ! Calculate factors that do depend on parameters, but not on RIJ
187:         CALL CALC_FCTS(X(NATOMS*3-8:NATOMS*3-6))137:         CALL CALC_FCTS(X(NATOMS*3-8:NATOMS*3-7))
188: 138: 
189:         IF (.NOT. PERIODIC) THEN139:         IF (.NOT. PERIODIC) THEN
190:             IF (LJ_GAUSS_MODE .EQ. 0) THEN140:             IF (LJ_GAUSS_MODE .EQ. 0) THEN
191:                 ! Normal cluster141:                 ! Normal cluster
192:                 NMOL = NATOMS142:                 NMOL = NATOMS
193:                 DO I = 1, NMOL-1 ! Outer loop over particles143:                 DO I = 1, NMOL-1 ! Outer loop over particles
194:                     DO J = I+1, NMOL ! Inner loop over particles144:                     DO J = I+1, NMOL ! Inner loop over particles
195:                         ! Get the particle-particle distance145:                         ! Get the particle-particle distance
196:                         RIJ(:) = X(3*I-2:3*I) - X(3*J-2:3*J)146:                         RIJ(:) = X(3*I-2:3*I) - X(3*J-2:3*J)
197:                         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))147:                         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))
250:                                           ENERGY, GTEST, &200:                                           ENERGY, GTEST, &
251:                                           GRAD(3*NATOMS-2:3*NATOMS))201:                                           GRAD(3*NATOMS-2:3*NATOMS))
252:                     END DO ! Inner loop over particles202:                     END DO ! Inner loop over particles
253:                 END DO ! Outer loop over particles203:                 END DO ! Outer loop over particles
254: 204: 
255:             CASE(2)205:             CASE(2)
256:                 ! Triclinic varying lattice206:                 ! Triclinic varying lattice
257:                 ! Reset all particles to within the box207:                 ! Reset all particles to within the box
258:                 CALL PERIODIC_RESET(X(:))208:                 CALL PERIODIC_RESET(X(:))
259: 209: 
260:                 ! Calculate box derivative factors 
261:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), & 
262:                                     X(3*NATOMS-2:3*NATOMS  )) 
263:  
264: !                WRITE(MYUNIT, *) 'Box lengths: ', X(3*NATOMS-2:3*NATOMS) 
265: !                WRITE(MYUNIT, *) 'Box angles: ', X(3*NATOMS-5:3*NATOMS-3) 
266: !                WRITE(MYUNIT, *) 'CELL_RANGE(:) ', BD%CELL_RANGE(:) 
267: !                WRITE(MYUNIT, *) 'V = ', BD%V 
268:  
269:                 ! Check whether the unit cell angles are physically possible.210:                 ! Check whether the unit cell angles are physically possible.
270:                 ! Reject if not.211:                 ! Reject if not.
271:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN212:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN
272:                     CALL REJECT(ENERGY, GRAD)213:                     CALL REJECT(ENERGY, GRAD)
273: !                    WRITE(MYUNIT, *) 'Reject for CHECK_ANGLES' 
274:                     RETURN214:                     RETURN
275:                 END IF215:                 END IF
276: 216: 
277:                 ! Reject steps where the cell range has got bigger than the 217:                 ! Calculate box derivative factors
278:                 ! cutoff. Such unit cells are highly distorted and there are 218:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &
279:                 ! probably equivalent ones of a more usual shape. Allowing them219:                                     X(3*NATOMS-2:3*NATOMS  ))
280:                 ! leads to very slow run times.220: 
281:                 IF (.NOT. ALL(BD%CELL_RANGE .LE. IMAGE_CUTOFF)) THEN221:                 ! Reject steps where the cell range has got bigger than 20. Such
 222:                 ! unit cells are highly distorted and there are probably
 223:                 ! equivalent ones of a more usual shape. Allowing them leads to
 224:                 ! very slow run times.
 225:                 IF (.NOT. ALL(BD%CELL_RANGE .LE. 20)) THEN
282:                     CALL REJECT(ENERGY, GRAD)226:                     CALL REJECT(ENERGY, GRAD)
283: !                    WRITE(MYUNIT, *) 'Reject for CELL_RANGE' 
284:                     RETURN227:                     RETURN
285:                 END IF228:                 END IF
286: 229: 
 230: !                WRITE(MYUNIT, *) 'Box lengths: ', X(3*NATOMS-2:3*NATOMS)
 231: !                WRITE(MYUNIT, *) 'Box angles: ', X(3*NATOMS-5:3*NATOMS-3)
 232: !                WRITE(MYUNIT, *) 'CELL_RANGE(:) ', BD%CELL_RANGE(:)
 233: 
287:                 NMOL = NATOMS-2234:                 NMOL = NATOMS-2
288:                 ! We need to include self-interactions of particles in different235:                 ! We need to include self-interactions of particles in different
289:                 ! unit cells236:                 ! unit cells
290:                 DO I = 1, NMOL! Outer loop over particles237:                 DO I = 1, NMOL! Outer loop over particles
291:                     DO J = I, NMOL ! Inner loop over particles238:                     DO J = I, NMOL ! Inner loop over particles
292:                         ! Add energy and gradients239:                         ! Add energy and gradients
293:                         CALL LJ_GAUSS_PER_TRI(I, J, X(:), GRAD(3*I-2:3*I), &240:                         CALL LJ_GAUSS_PER_TRI(I, J, X(:), GRAD(3*I-2:3*I), &
294:                             GRAD(3*J-2:3*J), ENERGY, GTEST, &241:                             GRAD(3*J-2:3*J), ENERGY, GTEST, &
295:                             GRAD(3*NATOMS-5:3*NATOMS), BD)242:                             GRAD(3*NATOMS-5:3*NATOMS), BD)
296:                     END DO ! Inner loop over particles243:                     END DO ! Inner loop over particles
297:                 END DO ! Outer loop over particles244:                 END DO ! Outer loop over particles
298:  
299:                 ! Constrain the box volume to be greater than 0 with a WCA-style 
300:                 ! repulsion 
301:                 CALL CONSTRAIN_VOLUME(ENERGY, GRAD(3*NATOMS-5:3*NATOMS-3), BD, & 
302:                                       GTEST) 
303:  
304:             CASE(3)245:             CASE(3)
305:                 ! Triclinic varying lattice, with varying parameters246:                 ! Triclinic varying lattice, with varying parameters
306:                 ! Reset all particles to within the box247:                 ! Reset all particles to within the box
307:                 CALL PERIODIC_RESET(X(:))248:                 CALL PERIODIC_RESET(X(:))
308: 249: 
309:                 ! Calculate box derivative factors 
310:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), & 
311:                                     X(3*NATOMS-2:3*NATOMS  )) 
312:  
313:                 ! Check whether the unit cell angles are physically possible.250:                 ! Check whether the unit cell angles are physically possible.
314:                 ! If we have a disallowed unit cell, we set the energy to very251:                 ! If we have a disallowed unit cell, we set the energy to very
315:                 ! high and the gradients to very small, so the step is252:                 ! high and the gradients to very small, so the step is
316:                 ! 'converged' but gets rejected.253:                 ! 'converged' but gets rejected.
317:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN254:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN
318:                     CALL REJECT(ENERGY, GRAD)255:                     ENERGY = 1.D20
 256:                     GRAD(:) = 1.D-20
 257:                     RETURN
319:                 END IF258:                 END IF
320: 259: 
321:                 ! Reject steps where the cell range has got bigger than the 260:                 ! Calculate box derivative factors
322:                 ! cutoff. Such unit cells are highly distorted and there are 261:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &
323:                 ! probably equivalent ones of a more usual shape. Allowing them262:                                     X(3*NATOMS-2:3*NATOMS  ))
324:                 ! leads to very slow run times.263: 
325:                 IF (.NOT. ALL(BD%CELL_RANGE .LE. IMAGE_CUTOFF)) THEN264:                 ! Reject steps where the cell range has got bigger than 20. Such
 265:                 ! unit cells are highly distorted and there are probably
 266:                 ! equivalent ones of a more usual shape. Allowing them leads to
 267:                 ! very slow run times.
 268:                 IF (.NOT. ALL(BD%CELL_RANGE .LE. 20)) THEN
326:                     CALL REJECT(ENERGY, GRAD)269:                     CALL REJECT(ENERGY, GRAD)
327:                     RETURN270:                     RETURN
328:                 END IF271:                 END IF
329: 272: 
330:                 NMOL = NATOMS-3273:                 NMOL = NATOMS-3
331:                 ! We need to include self-interactions of particles in different274:                 ! We need to include self-interactions of particles in different
332:                 ! unit cells275:                 ! unit cells
333:                 DO I = 1, NMOL! Outer loop over particles276:                 DO I = 1, NMOL! Outer loop over particles
334:                     DO J = I, NMOL ! Inner loop over particles277:                     DO J = I, NMOL ! Inner loop over particles
335:                         ! Add energy and gradients278:                         ! Add energy and gradients
336:                         CALL LJ_GAUSS_PER_TRI(I, J, X(:), GRAD(3*I-2:3*I), &279:                         CALL LJ_GAUSS_PER_TRI(I, J, X(:), GRAD(3*I-2:3*I), &
337:                             GRAD(3*J-2:3*J), ENERGY, GTEST, &280:                             GRAD(3*J-2:3*J), ENERGY, GTEST, &
338:                             GRAD(3*NATOMS-5:3*NATOMS), BD, &281:                             GRAD(3*NATOMS-5:3*NATOMS), BD, &
339:                             GRAD(3*NATOMS-8:3*NATOMS-6))282:                             GRAD(3*NATOMS-8:3*NATOMS-7))
340:                     END DO ! Inner loop over particles283:                     END DO ! Inner loop over particles
341:                 END DO ! Outer loop over particles284:                 END DO ! Outer loop over particles
342: 285: 
343:                 ! Constrain the box volume to be greater than 0 with a WCA-style 
344:                 ! repulsion 
345:                 CALL CONSTRAIN_VOLUME(ENERGY, GRAD(3*NATOMS-5:3*NATOMS-3), BD, & 
346:                                       GTEST) 
347:  
348:                 ! Constrain the potential parameters with a harmonic repulsion286:                 ! Constrain the potential parameters with a harmonic repulsion
349:                 CALL CONSTRAIN_PARAMETERS(ENERGY, GRAD(3*NATOMS-8:3*NATOMS-6), &287:                 CALL CONSTRAIN_PARAMETERS(ENERGY, GRAD(3*NATOMS-8:3*NATOMS-7), &
350:                                           GTEST)288:                                           GTEST)
351:             END SELECT ! Mode selections289:             END SELECT ! Mode selections
352:         END IF ! Periodic290:         END IF ! Periodic
353: 291: 
354:         ! Apply the overall scaling292:         ! Apply the overall scaling
355:         ENERGY = ENERGY * SCL_FCT293:         ENERGY = ENERGY * SCL_FCT
356:         GRAD(:) = GRAD(:) * SCL_FCT294:         GRAD(:) = GRAD(:) * SCL_FCT
357: 295: 
358:         ! The gradients wrt the potential parameters already include SCL_FCT296:         ! The gradients wrt the potential parameters already include SCL_FCT
359:         IF (LJ_GAUSS_MODE .EQ. 3) THEN297:         IF (LJ_GAUSS_MODE .EQ. 3) THEN
360:             GRAD(3*NATOMS-8:3*NATOMS-6) = GRAD(3*NATOMS-8:3*NATOMS-6) / SCL_FCT298:             GRAD(3*NATOMS-8:3*NATOMS-7) = GRAD(3*NATOMS-8:3*NATOMS-7) / SCL_FCT
361:         END IF299:         END IF
362: 300: 
363:     END SUBROUTINE LJ_GAUSS301:     END SUBROUTINE LJ_GAUSS
364: 302: 
365: !-------------------------------------------------------------------------------303: !-------------------------------------------------------------------------------
366: !304: !
367: ! Rejects a step by setting the energy very high and the gradients very low.305: ! Rejects a step by setting the energy very high and the gradients very low.
368: ! This tricks L-BFGS into thinking it's converged, but the high energy of the306: ! This tricks L-BFGS into thinking it's converged, but the high energy of the
369: ! quench will ensure it is rejected.307: ! quench will ensure it is rejected.
370: ! ENERGY: system energy308: ! ENERGY: system energy
392: ! distance between particles.330: ! distance between particles.
393: ! PARAMS: the potential parameters when those are to be adjusted, otherwise331: ! PARAMS: the potential parameters when those are to be adjusted, otherwise
394: !         ignored332: !         ignored
395: !333: !
396: !-------------------------------------------------------------------------------334: !-------------------------------------------------------------------------------
397: 335: 
398:     SUBROUTINE CALC_FCTS(PARAMS)336:     SUBROUTINE CALC_FCTS(PARAMS)
399: 337: 
400:         IMPLICIT NONE338:         IMPLICIT NONE
401: 339: 
402:         DOUBLE PRECISION, INTENT(INOUT) :: PARAMS(3)340:         DOUBLE PRECISION :: PARAMS(2)
403: 341: 
404:         ! Adjust the potential parameters, if the mode requires it342:         ! Adjust the potential parameters, if the mode requires it
405:         IF (LJ_GAUSS_MODE .EQ. 3) THEN343:         IF (LJ_GAUSS_MODE .EQ. 3) THEN
406:             ! For each parameter, if it is varying, set it from the coordinate.344:             LJ_GAUSS_EPS = PARAMS(1)
407:             ! If it's not varying, reset the coordinate, which TAKESTEP will345:             LJ_GAUSS_R0 = PARAMS(2)
408:             ! have changed. 
409:             IF (BTEST(LJ_GAUSS_PARAMS, 0)) THEN 
410:                 LJ_GAUSS_EPS = PARAMS(1) 
411:             ELSE 
412:                 PARAMS(1) = LJ_GAUSS_EPS 
413:             END IF 
414:             IF (BTEST(LJ_GAUSS_PARAMS, 1)) THEN 
415:                 LJ_GAUSS_R0 = PARAMS(2) 
416:             ELSE 
417:                 PARAMS(2) = LJ_GAUSS_R0 
418:             END IF 
419:             IF (BTEST(LJ_GAUSS_PARAMS, 2)) THEN 
420:                 LJ_GAUSS_SIGMASQ = PARAMS(3) 
421:             ELSE 
422:                 PARAMS(3) = LJ_GAUSS_SIGMASQ 
423:             END IF 
424:         END IF346:         END IF
425: 347: 
426:         ! Overall scale factor348:         ! Overall scale factor
427:         SCL_FCT = 1.D0 / (SQRT(2.D0 * PI * LJ_GAUSS_SIGMASQ) * LJ_GAUSS_EPS + &349:         SCL_FCT = 1.D0 / (SQRT(2.D0 * PI * LJ_GAUSS_SIGMASQ) * LJ_GAUSS_EPS + &
428:                           2.D0**(5.D0/6.D0) * 12.D0/55.D0)350:                           2.D0**(5.D0/6.D0) * 12.D0/55.D0)
429: 351: 
430:         ! Factors for the Gaussian352:         ! Factors for the Gaussian
431:         EXP_RCUT = LJ_GAUSS_EPS * EXP(- (LJ_GAUSS_RCUT - LJ_GAUSS_R0)**2 / &353:         EXP_RCUT = LJ_GAUSS_EPS * EXP(- (LJ_GAUSS_RCUT - LJ_GAUSS_R0)**2 / &
432:                                         (2.D0*LJ_GAUSS_SIGMASQ))354:                                         (2.D0*LJ_GAUSS_SIGMASQ))
433:         PREF_RCUT = (LJ_GAUSS_RCUT - LJ_GAUSS_R0) / LJ_GAUSS_SIGMASQ355:         PREF_RCUT = (LJ_GAUSS_RCUT - LJ_GAUSS_R0) / LJ_GAUSS_SIGMASQ
596: 518: 
597:         IMPLICIT NONE519:         IMPLICIT NONE
598: 520: 
599:         ! Check we have sensible value for the mode521:         ! Check we have sensible value for the mode
600:         IF (LJ_GAUSS_MODE .NE. 1 .AND. LJ_GAUSS_MODE .NE. 0 .AND. &522:         IF (LJ_GAUSS_MODE .NE. 1 .AND. LJ_GAUSS_MODE .NE. 0 .AND. &
601:             LJ_GAUSS_MODE .NE. 2 .AND. LJ_GAUSS_MODE .NE. 3) THEN523:             LJ_GAUSS_MODE .NE. 2 .AND. LJ_GAUSS_MODE .NE. 3) THEN
602:             WRITE (MYUNIT, *) 'LJ_GAUSS> ERROR: mode must be 0, 1, 2 or 3'524:             WRITE (MYUNIT, *) 'LJ_GAUSS> ERROR: mode must be 0, 1, 2 or 3'
603:             STOP525:             STOP
604:         END IF526:         END IF
605: 527: 
606:         ! Check we have sensible values for the params 
607:         IF (LJ_GAUSS_MODE .EQ. 3) THEN 
608:             IF (LJ_GAUSS_PARAMS .LT. 1 .OR. LJ_GAUSS_PARAMS .GT. 7) THEN 
609:                 WRITE (MYUNIT, *) 'LJ_GAUSS> ERROR: params must be between 1', & 
610:                                   ' and 7' 
611:                 STOP 
612:             END IF 
613:         END IF           
614:  
615:         ! Calculate some factors for LJ that don't depend on RIJ or the528:         ! Calculate some factors for LJ that don't depend on RIJ or the
616:         ! potential parameters529:         ! potential parameters
617:         RCUTINV6 = (1.D0/LJ_GAUSS_RCUT)**6530:         RCUTINV6 = (1.D0/LJ_GAUSS_RCUT)**6
618:         RCUTINV12 = RCUTINV6**2531:         RCUTINV12 = RCUTINV6**2
619: 532: 
620:     END SUBROUTINE LJ_GAUSS_INITIALISE533:     END SUBROUTINE LJ_GAUSS_INITIALISE
621: 534: 
622: !-------------------------------------------------------------------------------535: !-------------------------------------------------------------------------------
623: !536: !
624: ! Takes a step, perturbing the Cartesian coordinates, the unit cell parameters 
625: ! and the potential parameters, according to the specified mode. 
626: ! NP: the index in the main COORDS array to take the step from 
627: ! 
628: !------------------------------------------------------------------------------- 
629:  
630:     SUBROUTINE LJ_GAUSS_TAKESTEP(NP) 
631:  
632:         ! COORDS: full set of Markov chain coordinates 
633:         ! MYUNIT: file unit of the main GMIN output file 
634:         ! NATOMS: number of coordinates 
635:         ! PERIODIC: whetehr periodic boundary conditions are used 
636:         ! PERCOLATET: whether percolation is used to constrain particles 
637:         ! RADIUS: container radius 
638:         ! STEP: maximum step size for this move 
639:         ! TMOVE: specifies which step to do a translational move on 
640:         ! TWOD: whetehr the system is two dimensional 
641:         USE COMMONS, ONLY: COORDS, MYUNIT, NATOMS, PERIODIC, PERCOLATET, &  
642:                            RADIUS, STEP, TMOVE, TWOD 
643:          
644:         ! VEC_LEN: function, returns length of a 3D vector 
645:         ! VEC_RANDOM: function, returns a randomly orientated 3D unit vector 
646:         USE VEC3, ONLY: VEC_LEN, VEC_RANDOM 
647:  
648:         IMPLICIT NONE 
649:  
650:         INTEGER, INTENT(IN) :: NP 
651:  
652:         ! X: local copy of the coordinates 
653:         ! RAND_STEP: unit vector in a random step direction 
654:         ! STEP_SIZE: random fraction of the max step size to move 
655:         ! NE_ANGLES: new set of unit cell angles, to be checked 
656:         DOUBLE PRECISION :: X(3*NATOMS), RAND_STEP(3), STEP_SIZE, NEW_ANGLES(3) 
657:  
658:         ! I: iteration index 
659:         ! NMOL: actual number of atoms 
660:         INTEGER :: I, NMOL 
661:  
662:         ! Set the local copy of the coordinates 
663:         X(:) = COORDS(:, NP) 
664:  
665:         ! Set NMOL 
666:         NMOL = NATOMS - LJ_GAUSS_MODE 
667:  
668:         ! Random translational steps for each of the particles 
669:         IF (TMOVE(NP)) THEN 
670:             DO I = 1, NMOL 
671:                 ! Generate either a 2D or 3D random vector 
672:                 ! The magnitude will be uniformly distributed in a circle or  
673:                 ! sphere with radius given by the value of STEP for this step 
674:                 IF (TWOD) THEN 
675:                     DO 
676:                         CALL RANDOM_NUMBER(RAND_STEP(1)) 
677:                         CALL RANDOM_NUMBER(RAND_STEP(2)) 
678:                         RAND_STEP(3) = 0.D0 
679:                         IF(VEC_LEN(RAND_STEP) .LE. 1.D0) EXIT 
680:                     END DO 
681:                     RAND_STEP(:) = RAND_STEP(:) * STEP(NP) 
682:                 ELSE ! 3D 
683:                     RAND_STEP(:) = VEC_RANDOM() 
684:                     CALL RANDOM_NUMBER(STEP_SIZE) 
685:                     ! Sqrt for random volume distribution 
686:                     RAND_STEP(:) = RAND_STEP(:) * SQRT(STEP_SIZE * STEP(NP)) 
687:                 END IF 
688:                 ! Add the geenrated ste onto the coordinates 
689:                 X(I*3-2:I*3) = X(I*3-2:I*3) + RAND_STEP(:) 
690:             END DO 
691:         END IF 
692:  
693:         ! If not periodic and we're using a radius container, bring particles in 
694:         IF (.NOT. PERIODIC .AND. .NOT. PERCOLATET) THEN 
695:             DO I = 1, NMOL 
696:                 IF (VEC_LEN(X(I*3-2:I*3)) .GT. RADIUS) THEN 
697:                     WRITE(MYUNIT, *) 'LJ_GAUSS_TAKESTEP> ', & 
698:                         'coord outside container, bringing in' 
699:                     X(I*3-2:I*3) = X(I*3-2:I*3) - & 
700:                         SQRT(RADIUS) * NINT(X(I*3-2:I*3) / SQRT(RADIUS)) 
701:                 END IF 
702:             END DO 
703:         END IF 
704:  
705:         ! Box length steps 
706:         IF (LJ_GAUSS_MODE .GE. 1) THEN 
707:             X(3*NATOMS-2:3*NATOMS) = X(3*NATOMS-2:3*NATOMS) + & 
708:                                      VEC_RANDOM() * MAX_LENGTH_STEP 
709:         END IF 
710:  
711:         ! Box angle steps 
712:         DO 
713:             IF (LJ_GAUSS_MODE .GE. 2) THEN 
714:                 NEW_ANGLES(:) = X(3*NATOMS-5:3*NATOMS-3) + & 
715:                                 VEC_RANDOM() * MAX_ANGLE_STEP 
716:             END IF 
717:  
718:             ! See whether the unit cell we've generated is valid 
719:             IF (CHECK_ANGLES(NEW_ANGLES(:))) THEN 
720:                 X(3*NATOMS-5:3*NATOMS-3) =  NEW_ANGLES(:) 
721:                 EXIT 
722:             END IF 
723:         END DO 
724:  
725:         ! Parameter steps 
726:         IF (LJ_GAUSS_MODE .EQ. 3) THEN 
727:             IF (BTEST(LJ_GAUSS_PARAMS, 0)) THEN 
728:                 CALL RANDOM_NUMBER(STEP_SIZE) 
729:                 X(3*NATOMS-8) = X(3*NATOMS-8) + & 
730:                                 (STEP_SIZE - 0.5D0) * 2.D0 * MAX_EPS_STEP 
731:             END IF 
732:             IF (BTEST(LJ_GAUSS_PARAMS, 1)) THEN 
733:                 CALL RANDOM_NUMBER(STEP_SIZE) 
734:                 X(3*NATOMS-7) = X(3*NATOMS-7) + & 
735:                                 (STEP_SIZE - 0.5D0) * 2.D0 * MAX_R0_STEP 
736:             END IF 
737:             IF (BTEST(LJ_GAUSS_PARAMS, 2)) THEN 
738:                 CALL RANDOM_NUMBER(STEP_SIZE) 
739:                 X(3*NATOMS-6) = X(3*NATOMS-6) + & 
740:                                 (STEP_SIZE - 0.5D0) * 2.D0 * MAX_SIGMASQ_STEP 
741:             END IF             
742:         END IF 
743:  
744:         COORDS(:, NP) = X(:) 
745:  
746:     END SUBROUTINE LJ_GAUSS_TAKESTEP 
747:  
748: !------------------------------------------------------------------------------- 
749: ! 
750: ! Resets all coordinates to within the unit cell537: ! Resets all coordinates to within the unit cell
751: ! The unit cell runs from 0 to 1 in all directions. Coordinates are later scaled538: ! The unit cell runs from 0 to 1 in all directions. Coordinates are later scaled
752: ! by the box lengths.539: ! by the box lengths.
753: ! COORDS: coordinates of the particles540: ! COORDS: coordinates of the particles
754: !541: !
755: !-------------------------------------------------------------------------------542: !-------------------------------------------------------------------------------
756:     SUBROUTINE PERIODIC_RESET(COORDS)543:     SUBROUTINE PERIODIC_RESET(COORDS)
757: 544: 
758:         ! BOXLX, BOXLY, BOXLZ: dimensions of box545:         ! BOXLX, BOXLY, BOXLZ: dimensions of box
759:         ! NATOMS: number of particles546:         ! NATOMS: number of particles
784:             ! They will be checked later571:             ! They will be checked later
785:             COORDS(3*NATOMS-5:3*NATOMS-3) = COORDS(3*NATOMS-5:3*NATOMS-3) - &572:             COORDS(3*NATOMS-5:3*NATOMS-3) = COORDS(3*NATOMS-5:3*NATOMS-3) - &
786:                 2.D0*PI*FLOOR(COORDS(3*NATOMS-5:3*NATOMS-3)/(2.D0*PI))573:                 2.D0*PI*FLOOR(COORDS(3*NATOMS-5:3*NATOMS-3)/(2.D0*PI))
787:         CASE(3)574:         CASE(3)
788:            COORDS(1:3*NATOMS-9) = COORDS(1:3*NATOMS-9) - &575:            COORDS(1:3*NATOMS-9) = COORDS(1:3*NATOMS-9) - &
789:                                    FLOOR(COORDS(1:3*NATOMS-9))576:                                    FLOOR(COORDS(1:3*NATOMS-9))
790:             ! Make the lattice angles modulo 2*pi577:             ! Make the lattice angles modulo 2*pi
791:             ! They will be checked later578:             ! They will be checked later
792:             COORDS(3*NATOMS-5:3*NATOMS-3) = COORDS(3*NATOMS-5:3*NATOMS-3) - &579:             COORDS(3*NATOMS-5:3*NATOMS-3) = COORDS(3*NATOMS-5:3*NATOMS-3) - &
793:                 2.D0*PI*FLOOR(COORDS(3*NATOMS-5:3*NATOMS-3)/(2.D0*PI))580:                 2.D0*PI*FLOOR(COORDS(3*NATOMS-5:3*NATOMS-3)/(2.D0*PI))
794: 581:             ! Reset the unused parameter corodinate to zero
 582:             COORDS(3*NATOMS-6) = 0.D0
795:         END SELECT583:         END SELECT
796: 584: 
797:     END SUBROUTINE PERIODIC_RESET585:     END SUBROUTINE PERIODIC_RESET
798: 586: 
799: !-------------------------------------------------------------------------------587: !-------------------------------------------------------------------------------
800: !588: !
801: ! Finds the energy and gradients for the periodic system589: ! Finds the energy and gradients for the periodic system
802: ! IDI, IDJ: id numbers of the particles590: ! IDI, IDJ: id numbers of the particles
803: ! COORDS_I: coordinates of particle I591: ! COORDS_I: coordinates of particle I
804: ! COORDS_J: coordinates of particle J592: ! COORDS_J: coordinates of particle J
952:         ! Derivatives of the middle row wrt to lengths740:         ! Derivatives of the middle row wrt to lengths
953:         DERIV(2, 2, 5) = S(3)741:         DERIV(2, 2, 5) = S(3)
954:         DERIV(2, 3, 6) = (C(1) - C(2) * C(3)) / S(3)742:         DERIV(2, 3, 6) = (C(1) - C(2) * C(3)) / S(3)
955:         ! Derivatives of the bottom row wrt angles743:         ! 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)744:         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)745:         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)746:         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 lengths747:         ! Derivatives of the bottom row wrt lengths
960:         DERIV(3, 3, 6) = V / S(3)748:         DERIV(3, 3, 6) = V / S(3)
961: 749: 
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 vectors750:         ! Calculate the lengths of the reciprocal lattice vectors
968:         CALC_BOX_DERIV%RECIP(:) = S(:) / (SIZES(:) * V)751:         CALC_BOX_DERIV%RECIP(:) = S(:) / (SIZES(:) * V)
969: 752: 
970:         ! Calculate the number of cells to check in each direction753:         ! Calculate the number of cells to check in each direction
971:         CALC_BOX_DERIV%CELL_RANGE(:) = CEILING(LJ_GAUSS_RCUT * &754:         CALC_BOX_DERIV%CELL_RANGE(:) = CEILING(LJ_GAUSS_RCUT * &
972:             CALC_BOX_DERIV%RECIP(:))755:             CALC_BOX_DERIV%RECIP(:))
973: 756: 
974:         ! Set the output757:         ! Set the output
975:         CALC_BOX_DERIV%ORTHOG(:,:) = ORTHOG(:,:)758:         CALC_BOX_DERIV%ORTHOG(:,:) = ORTHOG(:,:)
976:         CALC_BOX_DERIV%DERIV(:,:,:) = DERIV(:,:,:)759:         CALC_BOX_DERIV%DERIV(:,:,:) = DERIV(:,:,:)
977:         CALC_BOX_DERIV%V = V 
978: 760: 
979:     END FUNCTION CALC_BOX_DERIV761:     END FUNCTION CALC_BOX_DERIV
980: 762: 
981: !-------------------------------------------------------------------------------763: !-------------------------------------------------------------------------------
982: !764: !
983: ! Checks whether the unit cell angles are allowed. It is possible to generate765: ! Checks whether the unit cell angles are allowed. Itis possible to generate
984: ! sets of angles which do not correspond to a physically possible unit cell.766: ! sets of angles which do not correspond to a physically possible unit cell.
985: ! This would manifest by an imaginary unit cell volume. See767: ! This would manifest by an imaginary unit cell volume. See
986: ! https://journals.iucr.org/a/issues/2011/01/00/au5114/au5114.pdf768: ! https://journals.iucr.org/a/issues/2011/01/00/au5114/au5114.pdf
987: ! 'On the allowed values for the triclinic unit-cell angles',769: ! 'On the allowed values for the triclinic unit-cell angles',
988: ! J. Foadi and G. Evans, Acta Cryst. (2011) A67, 93–95770: ! J. Foadi and G. Evans, Acta Cryst. (2011) A67, 93–95
989: ! ANGLES: the three unit cell angles771: ! ANGLES: the three unit cell angles
990: !772: !
991: !-------------------------------------------------------------------------------773: !-------------------------------------------------------------------------------
992: 774: 
993:     PURE LOGICAL FUNCTION CHECK_ANGLES(ANGLES)775:     PURE LOGICAL FUNCTION CHECK_ANGLES(ANGLES)
1036:         USE COMMONS, ONLY: NATOMS818:         USE COMMONS, ONLY: NATOMS
1037: 819: 
1038:         IMPLICIT NONE820:         IMPLICIT NONE
1039: 821: 
1040:         INTEGER, INTENT(IN) :: IDI, IDJ822:         INTEGER, INTENT(IN) :: IDI, IDJ
1041:         DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS)823:         DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS)
1042:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_I(3), GRAD_J(3), ENERGY824:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_I(3), GRAD_J(3), ENERGY
1043:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_BOX(6)825:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_BOX(6)
1044:         LOGICAL, INTENT(IN) :: GTEST826:         LOGICAL, INTENT(IN) :: GTEST
1045:         TYPE(BOX_DERIV), INTENT(IN) :: BD827:         TYPE(BOX_DERIV), INTENT(IN) :: BD
1046:         DOUBLE PRECISION, INTENT(INOUT), OPTIONAL :: GRAD_PARAMS(3)828:         DOUBLE PRECISION, INTENT(INOUT), OPTIONAL :: GRAD_PARAMS(2)
1047: 829: 
1048:         ! RJ: particle coordinates in crystallographic space830:         ! RJ: particle coordinates in crystallographic space
1049:         ! RIJ: particle-particle vector in crystallographic space831:         ! RIJ: particle-particle vector in crystallographic space
1050:         ! YIJ: particle-particle vector in orthogonal space832:         ! YIJ: particle-particle vector in orthogonal space
1051:         ! MODYIJ: particle-particle distance in orthogonal space833:         ! MODYIJ: particle-particle distance in orthogonal space
1052:         ! DVDY: derivative of pair potential wrt MODYIJ834:         ! DVDY: derivative of pair potential wrt MODYIJ
1053:         ! TEMP_GRAD: temporary gradient store835:         ! TEMP_GRAD: temporary gradient store
1054:         DOUBLE PRECISION :: RJ(3), RIJ(3), YIJ(3), MODYIJ836:         DOUBLE PRECISION :: RJ(3), RIJ(3), YIJ(3), MODYIJ
1055:         DOUBLE PRECISION :: DVDY, TEMP_GRAD(6)837:         DOUBLE PRECISION :: DVDY, TEMP_GRAD(6)
1056: 838: 
1136: ! Calculates the contribution to the derivatives of the potential wrt the918: ! Calculates the contribution to the derivatives of the potential wrt the
1137: ! potential parameters LJ_GAUSS_EPSILON and LJ_GAUSS_R0919: ! potential parameters LJ_GAUSS_EPSILON and LJ_GAUSS_R0
1138: ! MODYIJ: distance between the particles for this contribution920: ! MODYIJ: distance between the particles for this contribution
1139: !921: !
1140: !-------------------------------------------------------------------------------922: !-------------------------------------------------------------------------------
1141: 923: 
1142:     PURE FUNCTION CALC_GRAD_PARAMS(MODYIJ)924:     PURE FUNCTION CALC_GRAD_PARAMS(MODYIJ)
1143: 925: 
1144:         IMPLICIT NONE926:         IMPLICIT NONE
1145: 927: 
1146:         DOUBLE PRECISION :: CALC_GRAD_PARAMS(3)928:         DOUBLE PRECISION :: CALC_GRAD_PARAMS(2)
1147:         DOUBLE PRECISION, INTENT(IN) :: MODYIJ929:         DOUBLE PRECISION, INTENT(IN) :: MODYIJ
1148: 930: 
1149:         ! DVDE: derivative of pair potential wrt LJ_GAUSS_EPS931:         ! DVDE: derivative of pair potential wrt LJ_GAUSS_EPS
1150:         ! DSFDE: derivative of SCL_FCT wrt LJ_GAUSS_EPS932:         ! DSFDE: derivative of SCL_FCT wrt LJ_GAUSS_EPS
1151:         ! DVDR0: derivative of pair potential wrt LJ_GAUSS_R0933:         ! DVDR0: derivative of pair potential wrt LJ_GAUSS_R0
1152:         ! DVDS: derivative of pair potential wrt LJ_GAUSS_SIGMASQ 
1153:         ! DSFDS: derivative of SCL_FCT wrt LJ_GAUSS_SIGMASQ 
1154:         ! PREF_RIJ, EXP_RIJ: factors for the Guass potential934:         ! PREF_RIJ, EXP_RIJ: factors for the Guass potential
1155:         ! V: the pair potential energy935:         DOUBLE PRECISION :: DVDE, DSFDE, DVDR0, PREF_RIJ, EXP_RIJ
1156:         DOUBLE PRECISION :: DVDE, DSFDE, DVDR0, DVDS, DSFDS 
1157:         DOUBLE PRECISION :: PREF_RIJ, EXP_RIJ, V 
1158:  
1159:         ! Initialise output 
1160:         CALC_GRAD_PARAMS(:) = 0.D0 
1161: 936: 
1162:         ! Gauss factors937:         ! Gauss factors
1163:         PREF_RIJ = (MODYIJ - LJ_GAUSS_R0) / LJ_GAUSS_SIGMASQ938:         PREF_RIJ = (MODYIJ - LJ_GAUSS_R0) / LJ_GAUSS_SIGMASQ
1164:         EXP_RIJ = LJ_GAUSS_EPS * EXP(- (MODYIJ - LJ_GAUSS_R0)**2 / &939:         EXP_RIJ = LJ_GAUSS_EPS * EXP(- (MODYIJ - LJ_GAUSS_R0)**2 / &
1165:                                        (2.D0*LJ_GAUSS_SIGMASQ))940:                                        (2.D0*LJ_GAUSS_SIGMASQ))
1166: 941: 
1167:         ! Current energy942:         ! Derivative of the pair potential wrt LJ_GAUSS_EPS just requires
1168:         V = CALC_ENERGY(MODYIJ)943:         ! dividing by LJ_GAUSS_EPS
 944:         DVDE = CALC_ENERGY_GAUSS(MODYIJ) / LJ_GAUSS_EPS
 945: 
 946:         ! Derivative of the scale factor wrt LJ_GAUSS_EPS
 947:         DSFDE = - SQRT(2 * PI * LJ_GAUSS_SIGMASQ) * SCL_FCT**2
 948: 
 949:         ! Now we have everything we need for the LJ_GAUSS_EPS derivative
 950:         CALC_GRAD_PARAMS(1) = SCL_FCT * DVDE + DSFDE * CALC_ENERGY(MODYIJ)
 951: 
 952:         ! Next we look at the derivative wrt LJ_GAUSS_R0. This is simpler
 953:         ! because SCL_FCT does not depend on LJ_GAUSS_R0, but DVDR0 is more
 954:         ! complicated than DVDE
 955:         ! Gaussian derivative
 956:         DVDR0 = - PREF_RIJ * EXP_RIJ
 957:         ! Correction to make potential got to zero at RCUT
 958:         DVDR0 = DVDR0 + PREF_RCUT * EXP_RCUT
 959:         ! Correction to make the derivate go to zero at RCUT
 960:         DVDR0 = DVDR0 + (1.D0 / LJ_GAUSS_SIGMASQ - PREF_RCUT**2) * EXP_RCUT * &
 961:                         (MODYIJ - LJ_GAUSS_RCUT)
 962:         ! Correction to make the second derivative go to zero at RCUT
 963:         DVDR0 = DVDR0 - (3.D0 * PREF_RCUT/LJ_GAUSS_SIGMASQ - PREF_RCUT**3) * &
 964:                         EXP_RCUT * 0.5D0 * (MODYIJ - LJ_GAUSS_RCUT)**2
1169: 965: 
1170:         IF (BTEST(LJ_GAUSS_PARAMS, 0)) THEN966:         ! Now we can do the LJ_GAUSS_R0 derivative
1171:             ! Derivative of the pair potential wrt LJ_GAUSS_EPS just requires967:         CALC_GRAD_PARAMS(2) = SCL_FCT * DVDR0
1172:             ! dividing by LJ_GAUSS_EPS 
1173:             DVDE = CALC_ENERGY_GAUSS(MODYIJ) / LJ_GAUSS_EPS 
1174:  
1175:             ! Derivative of the scale factor wrt LJ_GAUSS_EPS 
1176:             DSFDE = - SQRT(2.D0 * PI * LJ_GAUSS_SIGMASQ) * SCL_FCT**2 
1177:  
1178:             ! Now we have everything we need for the LJ_GAUSS_EPS derivative 
1179:             CALC_GRAD_PARAMS(1) = SCL_FCT * DVDE + DSFDE * V 
1180:         END IF 
1181:  
1182:         IF (BTEST(LJ_GAUSS_PARAMS, 1)) THEN 
1183:             ! Next we look at the derivative wrt LJ_GAUSS_R0. This is simpler 
1184:             ! because SCL_FCT does not depend on LJ_GAUSS_R0, but DVDR0 is more 
1185:             ! complicated than DVDE 
1186:             ! Gaussian derivative 
1187:             DVDR0 = - PREF_RIJ * EXP_RIJ 
1188:             ! Correction to make potential got to zero at RCUT 
1189:             DVDR0 = DVDR0 + PREF_RCUT * EXP_RCUT 
1190:             ! Correction to make the derivate go to zero at RCUT 
1191:             DVDR0 = DVDR0 + (1.D0 / LJ_GAUSS_SIGMASQ - PREF_RCUT**2) * & 
1192:                              EXP_RCUT * (MODYIJ - LJ_GAUSS_RCUT) 
1193:             ! Correction to make the second derivative go to zero at RCUT 
1194:             DVDR0 = DVDR0 - (3.D0 * PREF_RCUT/LJ_GAUSS_SIGMASQ - PREF_RCUT**3)*& 
1195:                              EXP_RCUT * 0.5D0 * (MODYIJ - LJ_GAUSS_RCUT)**2 
1196:  
1197:             ! Now we can do the LJ_GAUSS_R0 derivative 
1198:             CALC_GRAD_PARAMS(2) = SCL_FCT * DVDR0 
1199:         END IF 
1200:  
1201:         IF (BTEST(LJ_GAUSS_PARAMS, 2)) THEN 
1202:             ! Finally we do the derivative wrt LJ_GAUSS_SIGMASQ. This one is 
1203:             ! really complicated... 
1204:             ! Gaussian derivative 
1205:             DVDS = - PREF_RIJ**2 * EXP_RIJ * 0.5D0 
1206:             ! Correction to make potential got to zero at RCUT 
1207:             DVDS = DVDS + PREF_RCUT**2 * EXP_RCUT * 0.5D0 
1208:             ! Correction to make the derivative go to zero at RCUT 
1209:             DVDS = DVDS + PREF_RCUT * EXP_RCUT * (MODYIJ - LJ_GAUSS_RCUT) * & 
1210:                    (1.D0 / LJ_GAUSS_SIGMASQ - PREF_RCUT**2 * 0.5D0) 
1211:             ! Correction to make the second derivative go to zero at RCUT 
1212:             DVDS = DVDS + EXP_RCUT * 0.5D0 * (MODYIJ - LJ_GAUSS_RCUT)**2 * & 
1213:                    (1.D0 / LJ_GAUSS_SIGMASQ**2 + PREF_RCUT**4 * 0.5D0 - & 
1214:                     5.D0 * PREF_RCUT**2 / (2.D0 * LJ_GAUSS_SIGMASQ)) 
1215:  
1216:             ! Derivative of SCL_FCT wrt LJ_GAUSS_SIGMASQ 
1217:             DSFDS = - SCL_FCT**2 * PI * LJ_GAUSS_EPS / & 
1218:                       SQRT(2.D0 * PI * LJ_GAUSS_SIGMASQ) 
1219:  
1220:             ! Now we have all the necessary factors 
1221:             CALC_GRAD_PARAMS(3) = SCL_FCT * DVDS + DSFDS * V 
1222:         END IF 
1223: 968: 
1224:     END FUNCTION CALC_GRAD_PARAMS969:     END FUNCTION CALC_GRAD_PARAMS
1225: 970: 
1226: !-------------------------------------------------------------------------------971: !-------------------------------------------------------------------------------
1227: !972: !
1228: ! Repel the unit cell volume from approaching 0, which repels the unit cell from 
1229: ! adopting physically impossible angle combinations. Uses a WCA potential 
1230: ! ENERGY: energy of the system 
1231: ! GRAD_ANGLES: gradients for the box angles 
1232: ! BD: information on the box and its derivatives 
1233: ! GTEST: whether to calculate gradients 
1234: ! 
1235: !------------------------------------------------------------------------------- 
1236:  
1237:     SUBROUTINE CONSTRAIN_VOLUME(ENERGY, GRAD_ANGLES, BD, GTEST) 
1238:  
1239:         IMPLICIT NONE 
1240:  
1241:         DOUBLE PRECISION, INTENT(INOUT) :: ENERGY, GRAD_ANGLES(3) 
1242:         TYPE(BOX_DERIV), INTENT(IN) :: BD 
1243:         LOGICAL, INTENT(IN) :: GTEST 
1244:  
1245: !        WRITE (*, *) 'BD%V = ', BD%V 
1246: !        WRITE (*, *) 'BD%V_DERIV(:) = ', BD%V_DERIV(:) 
1247:  
1248:         ! Add a purely repulsive (away from zero) WCA energy term 
1249:         IF (BD%V .LT. V_SIGMA**(1.D0/6.D0)) THEN 
1250:             ENERGY = ENERGY + 4.D0 * V_EPS * & 
1251:                      ((V_SIGMA / BD%V)**(12) - (V_SIGMA / BD%V)**(6)) + V_EPS 
1252:  
1253:             IF (GTEST) THEN 
1254:                 ! Add the gradients, if necessary 
1255:                 GRAD_ANGLES(:) = GRAD_ANGLES(:) + 24.D0 * V_EPS / V_SIGMA * & 
1256:                     ((V_SIGMA / BD%V)**(7) - 2.D0 * (V_SIGMA / BD%V)**(13)) * & 
1257:                     BD%V_DERIV(:) 
1258:             END IF 
1259:         END IF 
1260:  
1261:     END SUBROUTINE CONSTRAIN_VOLUME 
1262:  
1263: !------------------------------------------------------------------------------- 
1264: ! 
1265: ! Impose limits on the allowed range of the potential parameters LJ_GAUSS_R0 and973: ! Impose limits on the allowed range of the potential parameters LJ_GAUSS_R0 and
1266: ! LJ_GAUSS_EPSILON, by adding a harmonic repulsion outside the allowed range.974: ! LJ_GAUSS_EPSILON, by adding a harmonic repulsion outside the allowed range.
1267: ! We divide this energy contribution by SCL_FCT (it will be multiplied later) so975: ! We divide this energy contribution by SCL_FCT (it will be multiplied later) so
1268: ! the constraint is indepenednt of SCL_FCT.976: ! the constraint is indepenednt of SCL_FCT.
1269: ! ENERGY: energy of the system977: ! ENERGY: energy of the system
1270: ! GRAD_PARAMS: gradients for the potential parameters978: ! GRAD_PARAMS: gradients for the potential parameters
1271: ! GTEST: whether to calculate gradients979: ! GTEST: whether to calculate gradients
1272: !980: !
1273: !-------------------------------------------------------------------------------981: !-------------------------------------------------------------------------------
1274:     SUBROUTINE CONSTRAIN_PARAMETERS(ENERGY, GRAD_PARAMS, GTEST)982:     SUBROUTINE CONSTRAIN_PARAMETERS(ENERGY, GRAD_PARAMS, GTEST)
1275: 983: 
1276:         IMPLICIT NONE984:         IMPLICIT NONE
1277: 985: 
1278:         DOUBLE PRECISION, INTENT(INOUT) :: ENERGY, GRAD_PARAMS(3)986:         DOUBLE PRECISION, INTENT(INOUT) :: ENERGY, GRAD_PARAMS(2)
1279:         LOGICAL, INTENT(IN) :: GTEST987:         LOGICAL, INTENT(IN) :: GTEST
1280: 988: 
 989:         ! EPS_LOWER, EPS_UPPER: bounds on LJ_GAUSS_EPSILON
 990:         ! R0_LOWER, R0_UPPER: bounds on LJ_GAUSS_R0
 991:         ! REPULSION: harmonic force constant
 992:         DOUBLE PRECISION, PARAMETER :: EPS_LOWER = 0.D0, EPS_UPPER = 5.D0
 993:         DOUBLE PRECISION, PARAMETER :: R0_LOWER = 1.D0, R0_UPPER = 2.1D0
 994:         DOUBLE PRECISION, PARAMETER :: REPULSION = 1.D4
 995: 
1281:         ! EN_CON: constraint contribution to the energy996:         ! EN_CON: constraint contribution to the energy
1282:         DOUBLE PRECISION :: EN_CON997:         DOUBLE PRECISION :: EN_CON
1283: 998: 
1284:         ! Initialise output variable999:         ! Initialise output variable
1285:         EN_CON = 0.D01000:         EN_CON = 0.D0
1286: 1001: 
1287:         ! Epsilon energy contribution1002:         ! Epsilon energy contribution
1288:         IF (BTEST(LJ_GAUSS_PARAMS, 0)) THEN1003:         EN_CON = EN_CON + 0.5D0 * REPULSION * &
1289:             EN_CON = EN_CON + 0.5D0 * EPS_REPULSION * &1004:         (MERGE(LJ_GAUSS_EPS-EPS_UPPER, 0.D0, LJ_GAUSS_EPS .GT. EPS_UPPER)**2 + &
1290:                 (MERGE(LJ_GAUSS_EPS-EPS_UPPER, 0.D0, &1005:          MERGE(LJ_GAUSS_EPS-EPS_LOWER, 0.D0, LJ_GAUSS_EPS .LT. EPS_LOWER)**2)
1291:                        LJ_GAUSS_EPS .GT. EPS_UPPER)**2 + & 
1292:                  MERGE(LJ_GAUSS_EPS-EPS_LOWER, 0.D0, & 
1293:                        LJ_GAUSS_EPS .LT. EPS_LOWER)**2) 
1294:             IF (GTEST) THEN 
1295:             ! Epsilon gradient contribution 
1296:                 GRAD_PARAMS(1) = GRAD_PARAMS(1) + EPS_REPULSION * & 
1297:                     (MERGE(LJ_GAUSS_EPS-EPS_UPPER, 0.D0, & 
1298:                            LJ_GAUSS_EPS .GT. EPS_UPPER) + & 
1299:                      MERGE(LJ_GAUSS_EPS-EPS_LOWER, 0.D0, & 
1300:                            LJ_GAUSS_EPS .LT. EPS_LOWER)) 
1301:             END IF 
1302:         END IF 
1303: 1006: 
1304:         ! R0 energy contribution1007:         ! R0 energy contribution
1305:         IF (BTEST(LJ_GAUSS_PARAMS, 1)) THEN1008:         EN_CON = EN_CON + 0.5D0 * REPULSION * &
1306:             EN_CON = EN_CON + 0.5D0 * R0_REPULSION * &1009:             (MERGE(LJ_GAUSS_R0-R0_UPPER, 0.D0, LJ_GAUSS_R0 .GT. R0_UPPER)**2 + &
1307:                 (MERGE(LJ_GAUSS_R0 - R0_UPPER, 0.D0, &1010:              MERGE(LJ_GAUSS_R0-R0_LOWER, 0.D0, LJ_GAUSS_R0 .LT. R0_LOWER)**2)
1308:                        LJ_GAUSS_R0 .GT. R0_UPPER)**2 + & 
1309:                  MERGE(LJ_GAUSS_R0 - R0_LOWER, 0.D0, & 
1310:                        LJ_GAUSS_R0 .LT. R0_LOWER)**2) 
1311:             IF (GTEST) THEN 
1312:                 ! R0 gradient contribution 
1313:                 GRAD_PARAMS(2) = GRAD_PARAMS(2) + R0_REPULSION * & 
1314:                     (MERGE(LJ_GAUSS_R0 - R0_UPPER, 0.D0, & 
1315:                           LJ_GAUSS_R0 .GT. R0_UPPER) + & 
1316:                      MERGE(LJ_GAUSS_R0 - R0_LOWER, 0.D0, & 
1317:                            LJ_GAUSS_R0 .LT. R0_LOWER)) 
1318:             END IF 
1319:         END IF 
1320: 1011: 
1321:         ! SIGMASQ energy contribution1012:         IF (GTEST) THEN
1322:         IF (BTEST(LJ_GAUSS_PARAMS, 2)) THEN1013:             ! Epsilon gradient contribution
1323:             EN_CON = EN_CON + 0.5D0 * SIGMASQ_REPULSION * &1014:             GRAD_PARAMS(1) = GRAD_PARAMS(1) + REPULSION * &
1324:                 (MERGE(LJ_GAUSS_SIGMASQ - SIGMASQ_UPPER, 0.D0, &1015:              (MERGE(LJ_GAUSS_EPS-EPS_UPPER, 0.D0, LJ_GAUSS_EPS .GT. EPS_UPPER)+&
1325:                        LJ_GAUSS_SIGMASQ .GT. SIGMASQ_UPPER)**2 + &1016:               MERGE(LJ_GAUSS_EPS-EPS_LOWER, 0.D0, LJ_GAUSS_EPS .LT. EPS_LOWER))
1326:                  MERGE(LJ_GAUSS_SIGMASQ - SIGMASQ_LOWER, 0.D0, &1017:             ! R0 gradient contribution
1327:                        LJ_GAUSS_SIGMASQ .LT. SIGMASQ_LOWER)**2)1018:             GRAD_PARAMS(2) = GRAD_PARAMS(2) + REPULSION * &
1328:             IF (GTEST) THEN1019:                (MERGE(LJ_GAUSS_R0-R0_UPPER, 0.D0, LJ_GAUSS_R0 .GT. R0_UPPER) + &
1329:                 ! SIGMASQ gradient contribution1020:                 MERGE(LJ_GAUSS_R0-R0_LOWER, 0.D0, LJ_GAUSS_R0 .LT. R0_LOWER))
1330:                 GRAD_PARAMS(3) = GRAD_PARAMS(3) + SIGMASQ_REPULSION * & 
1331:                     (MERGE(LJ_GAUSS_SIGMASQ - SIGMASQ_UPPER, 0.D0, & 
1332:                            LJ_GAUSS_SIGMASQ .GT. SIGMASQ_UPPER) + & 
1333:                      MERGE(LJ_GAUSS_SIGMASQ - SIGMASQ_LOWER, 0.D0, & 
1334:                            LJ_GAUSS_SIGMASQ .LT. SIGMASQ_LOWER)) 
1335:             END IF 
1336:         END IF1021:         END IF
1337: 1022: 
1338:         ! Scale the contribtion and add it to the energy1023:         ! Scale the contribtion and add it to the energy
1339:         ENERGY = ENERGY + EN_CON / SCL_FCT1024:         ENERGY = ENERGY + EN_CON / SCL_FCT
1340: 1025: 
1341:     END SUBROUTINE CONSTRAIN_PARAMETERS1026:     END SUBROUTINE CONSTRAIN_PARAMETERS
1342: 1027: 
1343: !-------------------------------------------------------------------------------1028: !-------------------------------------------------------------------------------
1344: !1029: !
1345: ! Outputs the coordinates to the file ljgauss.xyz1030: ! Outputs the coordinates to the file ljgauss.xyz


r32371/mc.F 2017-04-25 16:30:10.298732784 +0100 r32370/mc.F 2017-04-25 16:30:10.970741476 +0100
 29:       USE MODAMBER9, ONLY : MDSTEPT,CISARRAY1,CHIARRAY1,NOCISTRANSDNA,NOCISTRANSRNA, 29:       USE MODAMBER9, ONLY : MDSTEPT,CISARRAY1,CHIARRAY1,NOCISTRANSDNA,NOCISTRANSRNA,
 30:      &                      SETCHIRAL,AMCHPMAX,DOLIGMOVE,LIGMOVEFREQ, AMCHNMAX, ROTAMERT, 30:      &                      SETCHIRAL,AMCHPMAX,DOLIGMOVE,LIGMOVEFREQ, AMCHNMAX, ROTAMERT,
 31:      &                      E_IGB, E_BOND, E_ANGLE, E_DIHEDRAL, 31:      &                      E_IGB, E_BOND, E_ANGLE, E_DIHEDRAL,
 32:      &                      E_VDW, E_14_VDW, E_ELEC, E_14_ELEC, IGB, MACROIONT 32:      &                      E_VDW, E_14_VDW, E_ELEC, E_14_ELEC, IGB, MACROIONT
 33:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_WRITE_RESTART, AMBER12_WRITE_PDB, POT_ENE_REC_C 33:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_WRITE_RESTART, AMBER12_WRITE_PDB, POT_ENE_REC_C
 34:       USE CHIRALITY, ONLY : INIT_CHIRAL, INIT_CIS_TRANS 34:       USE CHIRALITY, ONLY : INIT_CHIRAL, INIT_CIS_TRANS
 35:       USE porfuncs 35:       USE porfuncs
 36:       USE AMHGLOBALS, ONLY: NMRES,OMOVI,AVEP,NUMPRO,IRES 36:       USE AMHGLOBALS, ONLY: NMRES,OMOVI,AVEP,NUMPRO,IRES
 37:       USE AMH_INTERFACES, ONLY:E_WRITE 37:       USE AMH_INTERFACES, ONLY:E_WRITE
 38:       USE ROTAMER 38:       USE ROTAMER
 39:       USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS_TAKESTEP 
 40:  39: 
 41:       IMPLICIT NONE 40:       IMPLICIT NONE
 42: #ifdef MPI 41: #ifdef MPI
 43:       INCLUDE 'mpif.h' 42:       INCLUDE 'mpif.h'
 44:       INTEGER MPIERR 43:       INTEGER MPIERR
 45:       LOGICAL HITANY,MPI_FINT 44:       LOGICAL HITANY,MPI_FINT
 46: #endif 45: #endif
 47:  46: 
 48:       INTEGER J1, NSUCCESS(NPAR), NFAIL(NPAR), NFAILT(NPAR), NSUCCESST(NPAR), J2, NSTEPS, JP, REJSTREAK, 47:       INTEGER J1, NSUCCESS(NPAR), NFAIL(NPAR), NFAILT(NPAR), NSUCCESST(NPAR), J2, NSTEPS, JP, REJSTREAK,
 49:      1        UNT, ITERATIONS, NSUPERCOUNT, NQTOT, JACCPREV, NREN, NLAST, NSTEPREN, BRUN,QDONE,JBEST(NPAR),GCJBEST(NPAR), 48:      1        UNT, ITERATIONS, NSUPERCOUNT, NQTOT, JACCPREV, NREN, NLAST, NSTEPREN, BRUN,QDONE,JBEST(NPAR),GCJBEST(NPAR),
982: 981: 
983: !-------------------------------------------982: !-------------------------------------------
984: ! hk286 - generalised THOMSON983: ! hk286 - generalised THOMSON
985: !-----------------------------------------------984: !-----------------------------------------------
986:                ELSE IF (GTHOMSONT) THEN985:                ELSE IF (GTHOMSONT) THEN
987:                 CALL TAKESTEPGTHOMSON ()986:                 CALL TAKESTEPGTHOMSON ()
988: 987: 
989:                ELSE IF (GAYBERNEDCT) THEN988:                ELSE IF (GAYBERNEDCT) THEN
990:                 CALL TAKESTEPGB(JP)989:                 CALL TAKESTEPGB(JP)
991: 990: 
992: ! jwrm2> We need a special step taking routine for the LJGAUSS, to deal 
993: !        getting non-physical unit cell and potential parameters 
994:                ELSE IF (LJ_GAUSST) THEN 
995:                 CALL LJ_GAUSS_TAKESTEP(JP) 
996:  
997:                ELSE IF (PYT) THEN991:                ELSE IF (PYT) THEN
998: !               seed the random number generator with system time + MYNODE (for MPI runs)992: !               seed the random number generator with system time + MYNODE (for MPI runs)
999:                   IF(RANDOMSEEDT) THEN993:                   IF(RANDOMSEEDT) THEN
1000:                      CALL DATE_AND_TIME(datechar,timechar,zonechar,values)994:                      CALL DATE_AND_TIME(datechar,timechar,zonechar,values)
1001:                      itime1= values(6)*60 + values(7)995:                      itime1= values(6)*60 + values(7)
1002:                      CALL SDPRND(itime1+MYNODE)996:                      CALL SDPRND(itime1+MYNODE)
1003:                   END IF997:                   END IF
1004:                 call py_takestep(jp)998:                 call py_takestep(jp)
1005:                 if (rigidmdt) call rigidmd_nvt999:                 if (rigidmdt) call rigidmd_nvt
1006: 1000: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0