hdiff output

r33470/compress.f90 2017-11-13 18:30:12.993637599 +0000 r33469/compress.f90 2017-11-13 18:30:14.105652305 +0000
  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/compress.f90' in revision 33469
  2: !   Copyright (C) 1999-2006 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: !  Add energy and gradient correction terms for compression 
 21: ! 
 22: SUBROUTINE COMPRESS(X,V,ENERGY,GTEST) 
 23:     USE COMMONS 
 24:     USE GENRIGID, ONLY: NRIGIDBODY, RIGIDINIT 
 25:     USE OPP_MOD, ONLY: OPP_MODE 
 26:     IMPLICIT NONE 
 27:  
 28:     LOGICAL GTEST 
 29:     INTEGER J1, J3, NMOL 
 30:     DOUBLE PRECISION X(3*NATOMS), DIST, V(3*NATOMS), ENERGY, XMASS, YMASS, ZMASS 
 31:  
 32:     IF (PERIODIC) RETURN 
 33:  
 34: ! jwrm2> for rigid body angle axis system, only the first half of X is the body position 
 35: !        Note, as written this is a harmonic potential acting on the centre of the rigid 
 36: !        body, NOT on each site of the rigid body  
 37:     NMOL = NATOMS 
 38:     IF (RIGID) NMOL = NATOMS/2 
 39:     IF (RIGIDINIT) NMOL = NRIGIDBODY 
 40:  
 41: ! jwrm2> Need to remove alchemical and unit cell degrees of freedom. 
 42:     IF (OPPT) THEN 
 43:         SELECT CASE (OPP_MODE) 
 44:         CASE(1) 
 45:             NMOL = NATOMS - 2 
 46:         CASE(2) 
 47:             NMOL = NATOMS - 3 
 48:         CASE(3) 
 49:             NMOL = NATOMS - 1 
 50:         END SELECT 
 51:     END IF 
 52:  
 53: ! Find centre of mass 
 54:     XMASS=0.0D0 
 55:     YMASS=0.0D0 
 56:     ZMASS=0.0D0 
 57:     DO J1=1,NMOL 
 58:         XMASS=XMASS+X(3*(J1-1)+1) 
 59:         YMASS=YMASS+X(3*(J1-1)+2) 
 60:         ZMASS=ZMASS+X(3*(J1-1)+3) 
 61:     ENDDO 
 62:     XMASS=XMASS/NMOL 
 63:     YMASS=YMASS/NMOL 
 64:     ZMASS=ZMASS/NMOL 
 65:  
 66:     IF (DEBUG) WRITE (MYUNIT,*) 'compress> Compressing' 
 67:  
 68: ! Loop over all bodies 
 69:     DO J1=1, NMOL 
 70:         J3=3*J1 
 71:         DIST=(X(J3-2)-XMASS)**2+(X(J3-1)-YMASS)**2+(X(J3)-ZMASS)**2 
 72:  
 73: ! Add energy 
 74:         ENERGY=ENERGY+K_COMP*DIST/2.0D0 
 75:  
 76:         IF (GTEST) THEN 
 77: ! Add gradients 
 78:             V(J3-2)=V(J3-2)+K_COMP*(X(J3-2)-XMASS) 
 79:             V(J3-1)=V(J3-1)+K_COMP*(X(J3-1)-YMASS) 
 80:             V(J3)=  V(J3)  +K_COMP*(X(J3)  -ZMASS) 
 81:         ENDIF 
 82:     ENDDO 
 83:  
 84: END 


r33470/keywords.f 2017-11-13 18:30:13.221640611 +0000 r33469/keywords.f 2017-11-13 18:30:14.337655374 +0000
 38:       USE PORFUNCS 38:       USE PORFUNCS
 39:       USE MYGA_PARAMS 39:       USE MYGA_PARAMS
 40:       USE BGUPMOD 40:       USE BGUPMOD
 41:       USE GLJYMOD 41:       USE GLJYMOD
 42:       USE CHIRO_MODULE, ONLY: CHIRO_SIGMA, CHIRO_MU, CHIRO_GAMMA, CHIRO_L 42:       USE CHIRO_MODULE, ONLY: CHIRO_SIGMA, CHIRO_MU, CHIRO_GAMMA, CHIRO_L
 43:       USE CONVEX_POLYHEDRA_MODULE, ONLY: INITIALISE_POLYHEDRA, K_COMPRESS, K_OVERLAP 43:       USE CONVEX_POLYHEDRA_MODULE, ONLY: INITIALISE_POLYHEDRA, K_COMPRESS, K_OVERLAP
 44:       USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS_MODE, LJ_GAUSS_RCUT, LJ_GAUSS_EPS, 44:       USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS_MODE, LJ_GAUSS_RCUT, LJ_GAUSS_EPS,
 45:      &                        LJ_GAUSS_R0, LJ_GAUSS_SIGMASQ, LJ_GAUSS_PARAMS, 45:      &                        LJ_GAUSS_R0, LJ_GAUSS_SIGMASQ, LJ_GAUSS_PARAMS,
 46:      &                        LJ_GAUSS_INITIALISE 46:      &                        LJ_GAUSS_INITIALISE
 47:       USE OPP_MOD, ONLY: OPP_MODE, OPP_RCUT, OPP_K, OPP_PHI, OPP_PARAMS, 47:       USE OPP_MOD, ONLY: OPP_MODE, OPP_RCUT, OPP_K, OPP_PHI, OPP_PARAMS,
 48:      &                   OPP_INITIALISE, OPP_K_LOWER, OPP_K_UPPER, OPP_RON 48:      &                   OPP_INITIALISE, OPP_K_LOWER, OPP_K_UPPER
 49:       USE MBPOLMOD, ONLY: MBPOLINIT 49:       USE MBPOLMOD, ONLY: MBPOLINIT
 50:       USE SWMOD, ONLY: SWINIT, MWINIT 50:       USE SWMOD, ONLY: SWINIT, MWINIT
 51:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_GET_COORDS 51:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_GET_COORDS
 52: !     &                                 AMBER12_ATOMSS, AMBER12_SETUP, 52: !     &                                 AMBER12_ATOMSS, AMBER12_SETUP,
 53: !     &                                 AMBER12_RESIDUES, 53: !     &                                 AMBER12_RESIDUES,
 54: !     &                                 POPULATE_ATOM_DATA 54: !     &                                 POPULATE_ATOM_DATA
 55:       USE CHIRALITY, ONLY : CIS_TRANS_TOL 55:       USE CHIRALITY, ONLY : CIS_TRANS_TOL
 56:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR 56:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR
 57:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS, 57:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS,
 58:      &     PARSE_MLJ_PARAMS 58:      &     PARSE_MLJ_PARAMS
6593:           CALL READF(LJ_GAUSS_EPS)6593:           CALL READF(LJ_GAUSS_EPS)
6594:           CALL READF(LJ_GAUSS_R0)6594:           CALL READF(LJ_GAUSS_R0)
6595:           CALL READF(LJ_GAUSS_SIGMASQ)6595:           CALL READF(LJ_GAUSS_SIGMASQ)
6596:           IF (LJ_GAUSS_MODE .EQ. 3) CALL READI(LJ_GAUSS_PARAMS)6596:           IF (LJ_GAUSS_MODE .EQ. 3) CALL READI(LJ_GAUSS_PARAMS)
6597:           CALL LJ_GAUSS_INITIALISE()6597:           CALL LJ_GAUSS_INITIALISE()
6598: 6598: 
6599:       ELSE IF (WORD.EQ.'OPP') THEN6599:       ELSE IF (WORD.EQ.'OPP') THEN
6600:           OPPT = .TRUE.6600:           OPPT = .TRUE.
6601:           CALL READI(OPP_MODE)6601:           CALL READI(OPP_MODE)
6602:           CALL READF(OPP_RCUT)6602:           CALL READF(OPP_RCUT)
6603:           CALL READF(OPP_RON) 
6604:           CALL READF(OPP_K)6603:           CALL READF(OPP_K)
6605:           CALL READF(OPP_PHI)6604:           CALL READF(OPP_PHI)
6606:           IF (OPP_MODE .EQ. 2 .OR. OPP_MODE .EQ. 3) THEN6605:           IF (OPP_MODE .EQ. 2 .OR. OPP_MODE .EQ. 3) THEN
6607:               CALL READI(OPP_PARAMS)6606:               CALL READI(OPP_PARAMS)
6608:               IF (BTEST(OPP_PARAMS, 0)) THEN6607:               IF (BTEST(OPP_PARAMS, 0)) THEN
6609:                 CALL READF(OPP_K_LOWER)6608:                 CALL READF(OPP_K_LOWER)
6610:                 CALL READF(OPP_K_UPPER)6609:                 CALL READF(OPP_K_UPPER)
6611:               END IF6610:               END IF
6612:           END IF6611:           END IF
6613:           CALL OPP_INITIALISE()6612:           CALL OPP_INITIALISE()


r33470/opp.f90 2017-11-13 18:30:13.449643627 +0000 r33469/opp.f90 2017-11-13 18:30:14.561658338 +0000
 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: ! an oscillating pair potential (OPP). The possibility of varying 21: ! an oscillating pair potential (OPP). The possibility of varying
 22: ! and optimising the potential parameters is included. 22: ! and optimising the potential parameters is included.
 23: ! Also included is a framework for a varying triclinic unit cell. This could be 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. 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_RON, OPP_K, OPP_PHI, OPP_PARAMS 30: PUBLIC :: OPP_MODE, OPP_RCUT, OPP_K, OPP_PHI, OPP_PARAMS
 31: PUBLIC :: OPP_K_LOWER, OPP_K_UPPER 31: PUBLIC :: OPP_K_LOWER, OPP_K_UPPER
 32: ! Public subroutines 32: ! Public functions
 33: PUBLIC :: OPP, OPP_INITIALISE, OPP_TAKESTEP, OPP_VIEW 33: PUBLIC :: OPP, OPP_INITIALISE, VIEW_OPP, OPP_TAKESTEP
 34: ! Private parameters 34: ! Private parameters
 35: PRIVATE :: V_EPS, V_SIGMA, IMAGE_CUTOFF 35: PRIVATE :: V_EPS, V_SIGMA, FD_RCUT, IMAGE_CUTOFF!, SD_RCUT
 36: PRIVATE :: MAX_LENGTH_STEP, MAX_ANGLE_STEP, MAX_K_STEP, MAX_PHI_STEP 36: PRIVATE :: MAX_LENGTH_STEP, MAX_ANGLE_STEP, MAX_K_STEP, MAX_PHI_STEP
 37: PRIVATE :: K_REPULSION 37: PRIVATE :: K_REPULSION
 38: ! Private subroutines 
 39: PRIVATE :: REJECT, CALC_FCTS, PERIODIC_RESET, OPP_PER, OPP_PER_TRI, HESS_CONTRIB 
 40: PRIVATE :: HESS_CONTRIB_PER, CALC_HESS_PARAMS, CALC_HESS_PARAMS_PER 
 41: PRIVATE :: CONSTRAIN_VOLUME, CONSTRAIN_PARAMETERS 
 42: ! Private functions 38: ! Private functions
 43: PRIVATE :: CALC_ENERGY, CALC_GRAD, CALC_SEC, CALC_V, CALC_DVDR, CALC_D2VDR2 39: PRIVATE :: PERIODIC_RESET, OPP_PER, OPP_PER_TRI, CALC_ENERGY
 44: PRIVATE :: CALC_XPLOR, CALC_DXPLORDR, CALC_D2XPLORDR2, CHECK_ANGLES 40: PRIVATE :: CALC_GRAD, CALC_SEC, CALC_FCTS, CONSTRAIN_PARAMETERS
 45: PRIVATE :: CALC_GRAD_PARAMS, CALC_DVDK, CALC_DVDPHI, CALC_D2VDK2, CALC_D2VDPHI2 41: PRIVATE :: CHECK_ANGLES, REJECT, CONSTRAIN_VOLUME
 46: PRIVATE :: CALC_D2VDKDPHI, CALC_D2VDRDK, CALC_D2VDRDPHI 42: 
 47:  43: ! K: frequency parameter, controls the number of wells
 48: ! OPP_K: frequency parameter, controls the number of wells 44: ! PHI: phase parameter, controls the position of the first minimum
 49: ! OPP_PHI: phase parameter, controls the position of the first minimum 45: ! RCUT: cutoff distance at which the energy, gradients and Hessian of
 50: ! OPP_RCUT: cutoff distance at which the energy and gradients of the potential go  46: !       the potential go smoothly to zero.
 51: !       smoothly to zero. 47: DOUBLE PRECISION :: OPP_K, OPP_PHI, OPP_RCUT
 52: ! RON: distance at which the XPLOR smoothing begins 
 53: DOUBLE PRECISION :: OPP_K, OPP_PHI, OPP_RCUT, OPP_RON 
 54:  48: 
 55: ! OPP_MODE: 0 standard minimisation 49: ! OPP_MODE: 0 standard minimisation
 56: !           1 the final triplet of coordinates are the unit cell length 50: !           1 the final triplet of coordinates are the unit cell length
 57: !             parameters and the second last triplet are the unit cell 51: !             parameters and the second last triplet are the unit cell
 58: !             angles, giving a triclinic unit cell, which will be optimised 52: !             angles, giving a triclinic unit cell, which will be optimised
 59: !           2 as for 1, but the third last triplet are the potential 53: !           2 as for 1, but the third last triplet are the potential
 60: !             parameters OPP_K, OPP_PHI and 0 and 54: !             parameters OPP_K, OPP_PHI and 0 and
 61: !             the potential parameters will be optimised, depending on the 55: !             the potential parameters will be optimised, depending on the
 62: !             value of OPP_PARAMS 56: !             value of OPP_PARAMS
 63: !           3 for clusters, but with the last set of coordinates as the 57: !           3 for clusters, but with the last set of coordinates as the
 64: !             potential parameters 58: !             potential parameters
 65: ! OPP_PARAMS: Specifies which of the three potential parameters will be 59: ! OPP_PARAMS: Specifies which of the three potential parameters will be
 66: !             optimised. An integer between 1 and 3. 60: !             optimised. An integer between 1 and 3.
 67: !           1 OPP_K will be optimised 61: !           1 OPP_K will be optimised
 68: !           2 OPP_PHI will be optimised 62: !           2 OPP_PHI will be optimised
 69: !           3 both parameters will be optimised 63: !           3 both parameters will be optimised
 70: INTEGER :: OPP_MODE, OPP_PARAMS 64: INTEGER :: OPP_MODE, OPP_PARAMS
 71:  65: 
  66: ! V_RCUT: the unmodified potential evaluated at the cutoff
  67: ! FD_RCUT: the first derivative evaluated at the cutoff
  68: ! SD_RCUT: the second derivative evaluated at the cutoff
  69: DOUBLE PRECISION :: V_RCUT, FD_RCUT!, SD_RCUT
  70: 
 72: ! PI: the usual meaning 71: ! PI: the usual meaning
 73: DOUBLE PRECISION, PARAMETER :: PI = 4*ATAN(1.D0) 72: DOUBLE PRECISION, PARAMETER :: PI = 4*ATAN(1.D0)
 74:  73: 
 75: ! IMAGE_CUTOFF: if more periodic images than this in any direction fall 74: ! IMAGE_CUTOFF: if more periodic images than this in any direction fall
 76: !               within the cutoff, the evaluation will not be carried 75: !               within the cutoff, the evaluation will not be carried
 77: !               out and the step will be rejected. 76: !               out and the step will be rejected.
 78: DOUBLE PRECISION, PARAMETER :: IMAGE_CUTOFF = 30 77: DOUBLE PRECISION, PARAMETER :: IMAGE_CUTOFF = 30
 79:  78: 
 80: ! Parameters specifying the largest steps in OPP_TAKESTEP 79: ! Parameters specifying the largest steps in OPP_TAKESTEP
 81: ! MAX_LENGTH_STEP: largest allowed step in the length of a lattice 80: ! MAX_LENGTH_STEP: largest allowed step in the length of a lattice
133:         LOGICAL, INTENT(IN) :: GTEST, STEST132:         LOGICAL, INTENT(IN) :: GTEST, STEST
134: 133: 
135:         ! I, J: indices for looping over particles134:         ! I, J: indices for looping over particles
136:         ! NMOL: actual number of particles135:         ! NMOL: actual number of particles
137:         INTEGER :: I, J, NMOL136:         INTEGER :: I, J, NMOL
138: 137: 
139:         ! RIJ: Interparticle vector138:         ! RIJ: Interparticle vector
140:         ! MODRIJ: distance between a pair of particles139:         ! MODRIJ: distance between a pair of particles
141:         ! DVDR: derivative of pair potential wrt MODRIJ140:         ! DVDR: derivative of pair potential wrt MODRIJ
142:         ! D2VDR2: D(DVDR/MODRIJ)/DMODRIJ * 1/MODRIJ141:         ! D2VDR2: D(DVDR/MODRIJ)/DMODRIJ * 1/MODRIJ
143:         DOUBLE PRECISION :: RIJ(3), MODRIJ, DVDR142:         DOUBLE PRECISION :: RIJ(3), MODRIJ, DVDR, D2VDR2
144: 143: 
145:         ! Factors for triclinic lattice144:         ! Factors for triclinic lattice
146:         TYPE(BOX_DERIV) :: BD145:         TYPE(BOX_DERIV) :: BD
147: 146: 
148: !        WRITE (MYUNIT, *) 'OPP> beginning coords'147: !        WRITE (MYUNIT, *) 'OPP> beginning coords'
149: !        DO I = 1, NATOMS148: !        DO I = 1, NATOMS
150: !            WRITE (MYUNIT, *) X(3*I-2:3*I)149: !            WRITE (MYUNIT, *) X(3*I-2:3*I)
151: !        END DO150: !        END DO
152: 151: 
153:         ! Initialise output variables152:         ! Initialise output variables
154:         ENERGY = 0.D0153:         ENERGY = 0.D0
155:         GRAD(:) = 0.D0154:         GRAD(:) = 0.D0
156:         IF (STEST) HESS(:,:) = 0.D0155:         IF (STEST) HESS(:,:) = 0.D0
157: 156: 
158:         IF (.NOT. PERIODIC) THEN157:         IF (.NOT. PERIODIC) THEN
159:             IF (OPP_MODE .EQ. 0) THEN158:             IF (OPP_MODE .EQ. 0) THEN
160:                 ! Normal cluster159:                 ! Normal cluster
 160:                 CALL CALC_FCTS()
161:                 NMOL = NATOMS161:                 NMOL = NATOMS
162:                 DO I = 1, NMOL-1 ! Outer loop over particles162:                 DO I = 1, NMOL-1 ! Outer loop over particles
163:                     DO J = I+1, NMOL ! Inner loop over particles163:                     DO J = I+1, NMOL ! Inner loop over particles
164: 164: 
165:                         ! Get the particle-particle distance165:                         ! Get the particle-particle distance
166:                         RIJ(:) = X(3*I-2:3*I) - X(3*J-2:3*J)166:                         RIJ(:) = X(3*I-2:3*I) - X(3*J-2:3*J)
167:                         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))167:                         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))
168: 168: 
169:                         ! Check against cutoff169:                         ! Check against cutoff
170:                         IF (MODRIJ .LT. OPP_RCUT) THEN170:                         IF (MODRIJ .LT. OPP_RCUT) THEN
222:                                           GTEST, STEST)222:                                           GTEST, STEST)
223:             ELSE ! mode not equal to zero223:             ELSE ! mode not equal to zero
224:                 WRITE (MYUNIT, *) 'OPP> ERROR, PERIODIC must be ', &224:                 WRITE (MYUNIT, *) 'OPP> ERROR, PERIODIC must be ', &
225:                                   'specified if mode is not 0'225:                                   'specified if mode is not 0'
226:                 STOP226:                 STOP
227:             END IF227:             END IF
228:         ELSE ! PERIODIC228:         ELSE ! PERIODIC
229:             SELECT CASE (OPP_MODE)229:             SELECT CASE (OPP_MODE)
230:             CASE(0)230:             CASE(0)
231:                 ! Periodic, fixed box231:                 ! Periodic, fixed box
 232:                 CALL CALC_FCTS()
232:                 ! Reset all particles to within the box233:                 ! Reset all particles to within the box
233:                 CALL PERIODIC_RESET(X(:))234:                 CALL PERIODIC_RESET(X(:))
234: 235: 
235:                 NMOL = NATOMS236:                 NMOL = NATOMS
236:                 ! We need to include self-interactions of particles in different237:                 ! We need to include self-interactions of particles in different
237:                 ! unit cells238:                 ! unit cells
238:                 DO I = 1, NMOL ! Outer loop over particles239:                 DO I = 1, NMOL ! Outer loop over particles
239:                     DO J = 1, NMOL ! Inner loop over particles240:                     DO J = 1, NMOL ! Inner loop over particles
240:                         ! Add energy and gradients241:                         ! Add energy and gradients
241:                         CALL OPP_PER(I, J, X(3*I-2:3*I), X(3*J-2:3*J), &242:                         CALL OPP_PER(I, J, X(3*I-2:3*I), X(3*J-2:3*J), &
242:                                      GRAD(3*I-2:3*I), ENERGY, GTEST, STEST)243:                                      GRAD(3*I-2:3*I), ENERGY, GTEST, STEST)
243:                     END DO ! Inner loop over particles244:                     END DO ! Inner loop over particles
244:                 END DO ! Outer loop over particles245:                 END DO ! Outer loop over particles
245: 246: 
246:             CASE(1)247:             CASE(1)
247:                 ! Triclinic varying lattice248:                 ! Triclinic varying lattice
 249:                 CALL CALC_FCTS()
248:                 ! Reset all particles to within the box250:                 ! Reset all particles to within the box
249:                 CALL PERIODIC_RESET(X(:))251:                 CALL PERIODIC_RESET(X(:))
250: 252: 
251: !                DO I = 1, NATOMS253: !                DO I = 1, NATOMS
252: !                    WRITE (MYUNIT, *) X(3*I-2:3*I)254: !                    WRITE (MYUNIT, *) X(3*I-2:3*I)
253: !               END DO255: !               END DO
254: 256: 
255:                 ! Calculate box derivative factors257:                 ! Calculate box derivative factors
256:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &258:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &
257:                                     X(3*NATOMS-2:3*NATOMS  ), OPP_RCUT)259:                                     X(3*NATOMS-2:3*NATOMS  ), OPP_RCUT)
378: ! Some book keeping of the potential parameters to deal with different modes380: ! Some book keeping of the potential parameters to deal with different modes
379: ! PARAMS: the potential parameters when those are to be adjusted, otherwise381: ! PARAMS: the potential parameters when those are to be adjusted, otherwise
380: !         ignored382: !         ignored
381: !383: !
382: !-------------------------------------------------------------------------------384: !-------------------------------------------------------------------------------
383: 385: 
384:     SUBROUTINE CALC_FCTS(PARAMS)386:     SUBROUTINE CALC_FCTS(PARAMS)
385: 387: 
386:         IMPLICIT NONE388:         IMPLICIT NONE
387: 389: 
388:         DOUBLE PRECISION, INTENT(INOUT) :: PARAMS(3)390:         DOUBLE PRECISION, INTENT(INOUT), OPTIONAL :: PARAMS(3)
 391: 
 392:         ! S, C: factors which occur repeatedly
 393:         DOUBLE PRECISION :: S, C
389: 394: 
390:         ! Adjust the potential parameters, if the mode requires it395:         ! Adjust the potential parameters, if the mode requires it
391:         ! For each parameter, if it is varying, set it from the coordinate.396:         IF (PRESENT(PARAMS)) THEN
392:         ! If it's not varying, reset the coordinate, which TAKESTEP will397:             ! For each parameter, if it is varying, set it from the coordinate.
393:         ! have changed.398:             ! If it's not varying, reset the coordinate, which TAKESTEP will
394:         IF (BTEST(OPP_PARAMS, 0)) THEN399:             ! have changed.
395:             OPP_K = PARAMS(1)400:             IF (BTEST(OPP_PARAMS, 0)) THEN
396:         ELSE401:                 OPP_K = PARAMS(1)
397:             PARAMS(1) = OPP_K402:             ELSE
398:         END IF403:                 PARAMS(1) = OPP_K
399:         IF (BTEST(OPP_PARAMS, 1)) THEN404:             END IF
400:             OPP_PHI = PARAMS(2)405:             IF (BTEST(OPP_PARAMS, 1)) THEN
401:         ELSE406:                 OPP_PHI = PARAMS(2)
402:             PARAMS(2) = OPP_PHI407:             ELSE
 408:                 PARAMS(2) = OPP_PHI
 409:             END IF
 410:             PARAMS(3) = 0.D0
403:         END IF411:         END IF
404:         PARAMS(3) = 0.D0412: 
 413:         ! Calculate factors
 414:         S = SIN(OPP_K * (OPP_RCUT - 1) - OPP_PHI)
 415:         C = COS(OPP_K * (OPP_RCUT - 1) - OPP_PHI)
 416: 
 417:         ! Unmodifed potential evaluated at the cutoff
 418:         V_RCUT = OPP_RCUT**(-15) + OPP_RCUT**(-3) * C
 419:         ! First derivative evaluated at the cutoff
 420:         FD_RCUT = - 15.D0 * OPP_RCUT**(-16) - OPP_K * OPP_RCUT**(-3) * S - &
 421:                     3.D0 * OPP_RCUT**(-4) * C
 422:         ! Second derivative evaluated at the cutoff
 423: !        SD_RCUT = 240.D0 * OPP_RCUT**(-17) - OPP_K**2 * OPP_RCUT**(-3) * C + &
 424: !                  6.D0 * OPP_K * OPP_RCUT**(-4) * S + 12.D0 * OPP_RCUT**(-5) * C
405: 425: 
406:     END SUBROUTINE CALC_FCTS426:     END SUBROUTINE CALC_FCTS
407: 427: 
408: !-------------------------------------------------------------------------------428: !-------------------------------------------------------------------------------
409: !429: !
410: ! Calculates the energy contribution for a pair of particles.430: ! Calculates the energy contribution for a pair of particles.
411: ! MODRIJ: distance beteen the particle pair431: ! MODRIJ: distance beteen the particle pair
412: !432: !
413: !-------------------------------------------------------------------------------433: !-------------------------------------------------------------------------------
414:     PURE DOUBLE PRECISION FUNCTION CALC_ENERGY(MODRIJ)434:     PURE DOUBLE PRECISION FUNCTION CALC_ENERGY(MODRIJ)
415: 435: 
 436:         ! PERIODIC: whether peiodic boundary conditions are used
 437:         USE COMMONS, ONLY: PERIODIC
 438: 
416:         IMPLICIT NONE439:         IMPLICIT NONE
417: 440: 
418:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ441:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ
419: 442: 
420:         CALC_ENERGY = CALC_V(MODRIJ) * CALC_XPLOR(MODRIJ)443:         ! ENERGY: working variable for energy
 444:         DOUBLE PRECISION :: ENERGY
421: 445: 
422:     END FUNCTION CALC_ENERGY446:         ! Unmodified potential
 447:         ENERGY = MODRIJ**(-15) + MODRIJ**(-3) * &
 448:                  COS(OPP_K * (MODRIJ - 1) - OPP_PHI)
423: 449: 
424: !-------------------------------------------------------------------------------450:         ! Term to make the potential go to zero at OPP_RCUT
425: !451:         ENERGY = ENERGY - V_RCUT
426: ! Calculates the pair potential gradient for a pair of particles. 
427: ! MODRIJ: distance beteen the particle pair 
428: ! 
429: !------------------------------------------------------------------------------- 
430:     PURE DOUBLE PRECISION FUNCTION CALC_GRAD (MODRIJ) 
431: 452: 
432:         IMPLICIT NONE453:         ! Term to make the first derivative go to zero at OPP_RCUT
 454:         ENERGY = ENERGY - (MODRIJ - OPP_RCUT) * FD_RCUT
433: 455: 
434:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ456:         ! Term to make the second derivative go to zero at OPP_RCUT
 457: !        ENERGY = ENERGY - 0.5D0 * (MODRIJ - OPP_RCUT)**2 * SD_RCUT
435: 458: 
436:         CALC_GRAD = CALC_DVDR(MODRIJ) * CALC_XPLOR(MODRIJ) + &459:         ! Set the output
437:             CALC_V(MODRIJ) * CALC_DXPLORDR(MODRIJ)460:         CALC_ENERGY = ENERGY
438: 461: 
439:     END FUNCTION CALC_GRAD462:     END FUNCTION CALC_ENERGY
440: 463: 
441: !-------------------------------------------------------------------------------464: !-------------------------------------------------------------------------------
442: !465: !
443: ! Calculates the pair potential second derivative for a pair of particles.466: ! Calculates the pair potential gradient for a pair of particles.
444: ! MODRIJ: distance beteen the particle pair467: ! MODRIJ: distance beteen the particle pair
445: !468: !
446: !-------------------------------------------------------------------------------469: !-------------------------------------------------------------------------------
447:     PURE DOUBLE PRECISION FUNCTION CALC_SEC (MODRIJ)470:     PURE DOUBLE PRECISION FUNCTION CALC_GRAD (MODRIJ)
 471: 
 472:         ! PERIODIC: whether peiodic boundary conditions are used
 473:         USE COMMONS, ONLY: PERIODIC
448: 474: 
449:         IMPLICIT NONE475:         IMPLICIT NONE
450: 476: 
451:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ477:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ
452: 478: 
453:         CALC_SEC = CALC_D2VDR2(MODRIJ) * CALC_XPLOR(MODRIJ) + &479:         ! DVDR: temporary store for the result
454:             2.D0 * CALC_DVDR(MODRIJ) * CALC_DXPLORDR(MODRIJ) + &480:         DOUBLE PRECISION :: DVDR
455:             CALC_V(MODRIJ) * CALC_D2XPLORDR2(MODRIJ) 
456: 481: 
457:     END FUNCTION CALC_SEC482:         ! Unmodified gradient
 483:         DVDR = - 15.D0 * MODRIJ**(-16) - &
 484:                OPP_K * MODRIJ**(-3) * SIN(OPP_K * (MODRIJ - 1) - OPP_PHI) - & 
 485:                3.D0 * MODRIJ**(-4) * COS(OPP_K * (MODRIJ - 1) - OPP_PHI)
458: 486: 
459: !-------------------------------------------------------------------------------487:         ! Term to make the first derivative go to zero at OPP_RCUT
460: !488:         DVDR = DVDR - FD_RCUT
461: ! Calculates the potential unmodified by XPLOR 
462: ! MODRIJ: distance beteen the particle pair 
463: ! 
464: !------------------------------------------------------------------------------- 
465:     PURE DOUBLE PRECISION FUNCTION CALC_V (MODRIJ) 
466: 489: 
467:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ490:         ! Term to make the second derivative go to zero at OPP_RCUT
 491: !        DVDR = DVDR - (MODRIJ - OPP_RCUT) * SD_RCUT
468: 492: 
469:         CALC_V = MODRIJ**(-15) + MODRIJ**(-3) * &493:         CALC_GRAD = DVDR
470:             COS(OPP_K * (MODRIJ - 1) - OPP_PHI) 
471: 494: 
472:     END FUNCTION CALC_V495:     END FUNCTION CALC_GRAD
473: 496: 
474: !-------------------------------------------------------------------------------497: !-------------------------------------------------------------------------------
475: !498: !
476: ! Calculates the gradient of the potential unmodified by XPLOR499: ! Calculates the pair potential second derivative for a pair of particles.
477: ! MODRIJ: distance beteen the particle pair500: ! MODRIJ: distance beteen the particle pair
478: !501: !
479: !-------------------------------------------------------------------------------502: !-------------------------------------------------------------------------------
480:     PURE DOUBLE PRECISION FUNCTION CALC_DVDR (MODRIJ)503:     PURE DOUBLE PRECISION FUNCTION CALC_SEC (MODRIJ)
481:  
482:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
483:  
484:         CALC_DVDR = - 15.D0 * MODRIJ**(-16) - & 
485:             OPP_K * MODRIJ**(-3) * SIN(OPP_K * (MODRIJ - 1) - OPP_PHI) - &  
486:             3.D0 * MODRIJ**(-4) * COS(OPP_K * (MODRIJ - 1) - OPP_PHI) 
487:  
488:     END FUNCTION CALC_DVDR 
489: 504: 
490: !-------------------------------------------------------------------------------505:         IMPLICIT NONE
491: ! 
492: ! Calculates the second derivative of the potential unmodified by XPLOR 
493: ! MODRIJ: distance beteen the particle pair 
494: ! 
495: !------------------------------------------------------------------------------- 
496:     PURE DOUBLE PRECISION FUNCTION CALC_D2VDR2 (MODRIJ) 
497: 506: 
498:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ507:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ
499: 508: 
500:         CALC_D2VDR2 = 240.D0 * MODRIJ**(-17) - &509:         ! Unmodified gradient
 510:         CALC_SEC = 240.D0 * MODRIJ**(-17) - &
501:             OPP_K**2 * MODRIJ**(-3) * COS(OPP_K * (MODRIJ - 1) - OPP_PHI) + & 511:             OPP_K**2 * MODRIJ**(-3) * COS(OPP_K * (MODRIJ - 1) - OPP_PHI) + & 
502:             6.D0 * OPP_K * MODRIJ**(-4) * SIN(OPP_K * (MODRIJ - 1) - OPP_PHI) +&512:             6.D0 * OPP_K * MODRIJ**(-4) * SIN(OPP_K * (MODRIJ - 1) - OPP_PHI) +&
503:             12.D0 * MODRIJ**(-5) * COS(OPP_K * (MODRIJ - 1) - OPP_PHI)513:             12.D0 * MODRIJ**(-5) * COS(OPP_K * (MODRIJ - 1) - OPP_PHI)
504: 514: 
505:     END FUNCTION CALC_D2VDR2515:     END FUNCTION CALC_SEC
506:  
507: !------------------------------------------------------------------------------- 
508: ! 
509: ! Calculates the XPLOR smoothing potential 
510: ! MODRIJ: distance beteen the particle pair 
511: ! 
512: !------------------------------------------------------------------------------- 
513:     PURE DOUBLE PRECISION FUNCTION CALC_XPLOR (MODRIJ) 
514:  
515:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
516:  
517:         IF (MODRIJ .LT. OPP_RON) THEN 
518:             CALC_XPLOR = 1.D0 
519:         ELSE IF (MODRIJ .GT. OPP_RCUT) THEN 
520:             CALC_XPLOR = 0.D0 
521:         ELSE 
522:             CALC_XPLOR = (OPP_RCUT**2 - MODRIJ**2)**2 * & 
523:                 (OPP_RCUT**2 + 2.D0 * MODRIJ**2 - 3.D0 * OPP_RON**2) / & 
524:                 (OPP_RCUT**2 - OPP_RON**2)**3 
525:         END IF 
526:  
527:     END FUNCTION CALC_XPLOR 
528:  
529: !------------------------------------------------------------------------------- 
530: ! 
531: ! Calculates the derivative of the XPLOR smoothing potential 
532: ! MODRIJ: distance beteen the particle pair 
533: ! 
534: !------------------------------------------------------------------------------- 
535:     PURE DOUBLE PRECISION FUNCTION CALC_DXPLORDR (MODRIJ) 
536:  
537:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
538:  
539:         IF (MODRIJ .LT. OPP_RON .OR. MODRIJ .GT. OPP_RCUT) THEN 
540:             CALC_DXPLORDR = 0.D0 
541:         ELSE 
542:             CALC_DXPLORDR = 12.D0 * MODRIJ * (OPP_RCUT**2 - MODRIJ**2) * & 
543:                 (OPP_RON**2 - MODRIJ**2) / (OPP_RCUT**2 - OPP_RON**2)**3 
544:         END IF 
545:  
546:     END FUNCTION CALC_DXPLORDR 
547:  
548: !------------------------------------------------------------------------------- 
549: ! 
550: ! Calculates the second derivative of the XPLOR smoothing potential 
551: ! MODRIJ: distance beteen the particle pair 
552: ! 
553: !------------------------------------------------------------------------------- 
554:     PURE DOUBLE PRECISION FUNCTION CALC_D2XPLORDR2 (MODRIJ) 
555:  
556:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
557:  
558:         IF (MODRIJ .LT. OPP_RON .OR. MODRIJ .GT. OPP_RCUT) THEN 
559:             CALC_D2XPLORDR2 = 0.D0 
560:         ELSE 
561:             CALC_D2XPLORDR2 = 12.D0 * & 
562:                 ((OPP_RCUT**2 - MODRIJ**2) * (OPP_RON**2 - MODRIJ**2) - & 
563:                 2.D0 * MODRIJ**2 * (OPP_RCUT**2 - MODRIJ**2) - & 
564:                 2.D0 * MODRIJ**2 * (OPP_RON**2 - MODRIJ**2)) / & 
565:                 (OPP_RCUT**2 - OPP_RON**2)**3 
566:         END IF 
567:  
568:     END FUNCTION CALC_D2XPLORDR2 
569: 516: 
570: !-------------------------------------------------------------------------------517: !-------------------------------------------------------------------------------
571: !518: !
572: ! Initialisation. Check values for OPP_MODE and OPP_PARAMS.519: ! Initialisation. Check values for OPP_MODE and OPP_PARAMS.
573: !520: !
574: !-------------------------------------------------------------------------------521: !-------------------------------------------------------------------------------
575:     SUBROUTINE OPP_INITIALISE()522:     SUBROUTINE OPP_INITIALISE()
576: 523: 
577:         ! MYUNIT: file unit for main GMIN output524:         ! MYUNIT: file unit for main GMIN output
578:         USE COMMONS, ONLY: MYUNIT525:         USE COMMONS, ONLY: MYUNIT
584:             WRITE (MYUNIT, *) 'OPP> ERROR: mode must be between 0 and 3'531:             WRITE (MYUNIT, *) 'OPP> ERROR: mode must be between 0 and 3'
585:             STOP532:             STOP
586:         END IF533:         END IF
587: 534: 
588:         ! Check we have sensible values for the params535:         ! Check we have sensible values for the params
589:         IF (OPP_MODE .EQ. 2 .OR. OPP_MODE .EQ. 3) THEN536:         IF (OPP_MODE .EQ. 2 .OR. OPP_MODE .EQ. 3) THEN
590:             IF (OPP_PARAMS .LT. 1 .OR. OPP_PARAMS .GT. 3) THEN537:             IF (OPP_PARAMS .LT. 1 .OR. OPP_PARAMS .GT. 3) THEN
591:                 WRITE (MYUNIT, *) 'OPP> ERROR: params must be between 1 and 3'538:                 WRITE (MYUNIT, *) 'OPP> ERROR: params must be between 1 and 3'
592:                 STOP539:                 STOP
593:             END IF540:             END IF
594:         END IF541:         END IF          
595:  
596:         ! Check the XPLOR cut-off values 
597:         IF (OPP_RCUT .LE. OPP_RON) THEN 
598:             WRITE (MYUNIT, *) 'OPP> ERROR: XPLOR switch on must be less than', & 
599:                               ' the cut-off' 
600:             STOP 
601:         END IF     
602: 542: 
603:     END SUBROUTINE OPP_INITIALISE543:     END SUBROUTINE OPP_INITIALISE
604: 544: 
605: !-------------------------------------------------------------------------------545: !-------------------------------------------------------------------------------
606: !546: !
607: ! Takes a step, perturbing the Cartesian coordinates, the unit cell parameters547: ! Takes a step, perturbing the Cartesian coordinates, the unit cell parameters
608: ! and the potential parameters, according to the specified mode.548: ! and the potential parameters, according to the specified mode.
609: ! NP: the index in the main COORDS array to take the step from549: ! NP: the index in the main COORDS array to take the step from
610: !550: !
611: !-------------------------------------------------------------------------------551: !-------------------------------------------------------------------------------
1110:                                 MATMUL(TRANSPOSE(BD%ORTHOG(:,:)), YIJ(:))1050:                                 MATMUL(TRANSPOSE(BD%ORTHOG(:,:)), YIJ(:))
1111:                             GRAD_I(:) = GRAD_I(:) + TEMP_GRAD(1:3)1051:                             GRAD_I(:) = GRAD_I(:) + TEMP_GRAD(1:3)
1112:                             GRAD_J(:) = GRAD_J(:) - TEMP_GRAD(1:3)1052:                             GRAD_J(:) = GRAD_J(:) - TEMP_GRAD(1:3)
1113: 1053: 
1114:                             ! Lattice parameter derivatives1054:                             ! Lattice parameter derivatives
1115:                             DO I = 1, 6 ! Loop over the different box parameters1055:                             DO I = 1, 6 ! Loop over the different box parameters
1116:                                 TEMP = DVDY * &1056:                                 TEMP = DVDY * &
1117:                                     DOT_PRODUCT(YIJ(:), &1057:                                     DOT_PRODUCT(YIJ(:), &
1118:                                     MATMUL(BD%DERIV(:,:,I), RIJ(:)))1058:                                     MATMUL(BD%DERIV(:,:,I), RIJ(:)))
1119:                                 GRAD_BOX(I) = GRAD_BOX(I) + TEMP1059:                                 GRAD_BOX(I) = GRAD_BOX(I) + TEMP
 1060:                                 IF (I .EQ. 1 .AND. .NOT. STEST) THEN
 1061: !                                    WRITE (*, *) 'IDI = ', IDI, 'IDJ = ', IDJ, 'Gradient Contribution = ', TEMP
 1062:                                END IF
1120:                             END DO1063:                             END DO
1121: 1064: 
1122:                             ! If necessary, calculate the parameter derivatives1065:                             ! If necessary, calculate the parameter derivatives
1123:                             ! Divide gradient by 2 for images of the same atom,1066:                             ! Divide gradient by 2 for images of the same atom,
1124:                             ! to avoid double counting1067:                             ! to avoid double counting
1125:                             IF(PRESENT(GRAD_PARAMS) .AND. OPP_MODE .EQ. 2)&1068:                             IF(PRESENT(GRAD_PARAMS) .AND. OPP_MODE .EQ. 2)&
1126:                                 GRAD_PARAMS(:) = GRAD_PARAMS(:) + &1069:                                 GRAD_PARAMS(:) = GRAD_PARAMS(:) + &
1127:                                     MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * &1070:                                     MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * &
1128:                                     CALC_GRAD_PARAMS(MODYIJ)1071:                                     CALC_GRAD_PARAMS(MODYIJ)
1129: 1072: 
1142: !-------------------------------------------------------------------------------1085: !-------------------------------------------------------------------------------
1143: !1086: !
1144: ! Calculates the second derivative contributions for a cluster1087: ! Calculates the second derivative contributions for a cluster
1145: ! RIJ: vector between particles in fractional coordinates1088: ! RIJ: vector between particles in fractional coordinates
1146: ! BD: box derivative information1089: ! BD: box derivative information
1147: ! IDI, IDJ: indices of the atom pair1090: ! IDI, IDJ: indices of the atom pair
1148: !1091: !
1149: !-------------------------------------------------------------------------------1092: !-------------------------------------------------------------------------------
1150:     SUBROUTINE HESS_CONTRIB(RIJ, IDI, IDJ)1093:     SUBROUTINE HESS_CONTRIB(RIJ, IDI, IDJ)
1151: 1094: 
 1095:         ! Number of atoms, including unit cell parameters
 1096:         USE COMMONS, ONLY: NATOMS
 1097: 
1152:         ! HESS: the Hessian matrix1098:         ! HESS: the Hessian matrix
1153:         USE MODHESS, ONLY: HESS1099:         USE MODHESS, ONLY: HESS
1154: 1100: 
1155:         IMPLICIT NONE1101:         IMPLICIT NONE
1156: 1102: 
1157:         DOUBLE PRECISION, INTENT(IN) :: RIJ(3)1103:         DOUBLE PRECISION, INTENT(IN) :: RIJ(3)
1158:         INTEGER, INTENT(IN) :: IDI, IDJ1104:         INTEGER, INTENT(IN) :: IDI, IDJ
1159: 1105: 
1160:         ! MODRIJ: particle-particle distance in Cartesian space1106:         ! MODRIJ: particle-particle distance in Cartesian space
1161:         ! DVDR: derivative of pair potential wrt MODYIJ1107:         ! DVDR: derivative of pair potential wrt MODYIJ
1162:         ! D2VDR2: D(DVDR/MODRIJ)/DMODRIJ * 1/MODRIJ1108:         ! D2VDR2: D(DVDR/MODRIJ)/DMODRIJ * 1/MODRIJ
1163:         DOUBLE PRECISION :: MODRIJ, DVDR, D2VDR21109:         ! TEMP, TEMP_GRAD: temporary calculation store
 1110:         DOUBLE PRECISION :: MODRIJ, DVDR, D2VDR2, TEMP, TEMP_GRAD(3)
1164: 1111: 
1165:         ! I, J: loop indices1112:         ! I, J: loop indices
1166:         INTEGER :: I, J1113:         INTEGER :: I, J
1167: 1114: 
1168:         ! Calculate some factors1115:         ! Calculate some factors
1169:         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))1116:         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))
1170:         DVDR = CALC_GRAD(MODRIJ) / MODRIJ1117:         DVDR = CALC_GRAD(MODRIJ) / MODRIJ
1171:         D2VDR2 = (CALC_SEC(MODRIJ) - DVDR) / MODRIJ**21118:         D2VDR2 = (CALC_SEC(MODRIJ) - DVDR) / MODRIJ**2
1172: 1119: 
1173:         DO I = 1, 3 ! outer loop over x, y, z1120:         DO I = 1, 3 ! outer loop over x, y, z
1278:         END DO ! end loop over unit cell parameters1225:         END DO ! end loop over unit cell parameters
1279: 1226: 
1280:         ! Pure unit cell parameters1227:         ! Pure unit cell parameters
1281:         DO I = 1, 6 ! outer loop over cell parameters1228:         DO I = 1, 6 ! outer loop over cell parameters
1282:             ! Diagonals1229:             ! Diagonals
1283:             TEMP = D2VDY2 * &1230:             TEMP = D2VDY2 * &
1284:                 DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,I), RIJ(:)))**2 + &1231:                 DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,I), RIJ(:)))**2 + &
1285:                 DVDY * (DOT_PRODUCT(MATMUL(BD%DERIV(:,:,I), RIJ(:)), &1232:                 DVDY * (DOT_PRODUCT(MATMUL(BD%DERIV(:,:,I), RIJ(:)), &
1286:                 MATMUL(BD%DERIV(:,:,I), RIJ(:))) + &1233:                 MATMUL(BD%DERIV(:,:,I), RIJ(:))) + &
1287:                 DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV2(:,:,I,I), RIJ(:))))1234:                 DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV2(:,:,I,I), RIJ(:))))
1288:             HESS(3*NATOMS-6+I,3*NATOMS-6+I) = HESS(3*NATOMS-6+I,3*NATOMS-6+I) +&1235:             HESS(3*NATOMS-6+I,3*NATOMS-6+I) = HESS(3*NATOMS-6+I,3*NATOMS-6+I) + TEMP
1289:                 TEMP 
1290:             DO J = I + 1, 6 ! inner loop over cell parameters1236:             DO J = I + 1, 6 ! inner loop over cell parameters
1291:                 TEMP = D2VDY2 * &1237:                 TEMP = D2VDY2 * &
1292:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,I), RIJ(:))) * &1238:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,I), RIJ(:))) * &
1293:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,J), RIJ(:))) + &1239:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,J), RIJ(:))) + &
1294:                     DVDY * (DOT_PRODUCT(MATMUL(BD%DERIV(:,:,I), RIJ(:)), &1240:                     DVDY * (DOT_PRODUCT(MATMUL(BD%DERIV(:,:,I), RIJ(:)), &
1295:                     MATMUL(BD%DERIV(:,:,J), RIJ(:))) + &1241:                     MATMUL(BD%DERIV(:,:,J), RIJ(:))) + &
1296:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV2(:,:,I,J), RIJ(:))))1242:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV2(:,:,I,J), RIJ(:))))
1297:                 HESS(3*NATOMS-6+I,3*NATOMS-6+J) = &1243:                 HESS(3*NATOMS-6+I,3*NATOMS-6+J) = &
1298:                     HESS(3*NATOMS-6+I,3*NATOMS-6+J) + TEMP1244:                     HESS(3*NATOMS-6+I,3*NATOMS-6+J) + TEMP
1299:                 HESS(3*NATOMS-6+J,3*NATOMS-6+I) = &1245:                 HESS(3*NATOMS-6+J,3*NATOMS-6+I) = &
1311: !-------------------------------------------------------------------------------1257: !-------------------------------------------------------------------------------
1312: !1258: !
1313: ! Calculates the contribution to the derivatives of the potential wrt the1259: ! Calculates the contribution to the derivatives of the potential wrt the
1314: ! potential parameters OPP_K and OPP_PHI1260: ! potential parameters OPP_K and OPP_PHI
1315: ! MODRIJ: distance between the particles for this contribution1261: ! MODRIJ: distance between the particles for this contribution
1316: !1262: !
1317: !-------------------------------------------------------------------------------1263: !-------------------------------------------------------------------------------
1318: 1264: 
1319:     PURE FUNCTION CALC_GRAD_PARAMS(MODRIJ)1265:     PURE FUNCTION CALC_GRAD_PARAMS(MODRIJ)
1320: 1266: 
 1267:         USE COMMONS, ONLY: PERIODIC
 1268: 
1321:         IMPLICIT NONE1269:         IMPLICIT NONE
1322: 1270: 
1323:         DOUBLE PRECISION :: CALC_GRAD_PARAMS(3)1271:         DOUBLE PRECISION :: CALC_GRAD_PARAMS(3)
1324:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ1272:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ
1325: 1273: 
 1274:         ! S, C: sin and cos factor that appears repeatedly
 1275:         ! T: working store for the output
 1276:         DOUBLE PRECISION :: S, C
 1277:         DOUBLE PRECISION :: T(3)
 1278: 
1326:         ! Initialise output1279:         ! Initialise output
1327:         CALC_GRAD_PARAMS(:) = 0.D01280:         T(:) = 0.D0
 1281: 
 1282:         ! Calculate factors
 1283:         S = SIN(OPP_K * (OPP_RCUT - 1) - OPP_PHI)
 1284:         C = COS(OPP_K * (OPP_RCUT - 1) - OPP_PHI)
1328: 1285: 
1329:         IF (BTEST(OPP_PARAMS, 0)) THEN1286:         IF (BTEST(OPP_PARAMS, 0)) THEN
1330:             CALC_GRAD_PARAMS(1) = CALC_DVDK(MODRIJ) * CALC_XPLOR(MODRIJ)1287:             ! OPP_K derivative
 1288:             ! Unmodified potential term
 1289:             T(1) = - SIN(OPP_K * (MODRIJ-1) - OPP_PHI) * (MODRIJ-1) / MODRIJ**3
 1290:             ! Term to make the potential go to zero at OPP_RCUT
 1291:             T(1) = T(1) + S * (OPP_RCUT - 1) / OPP_RCUT **3
 1292:             ! Terms to make the first derivative go to zero at OPP_RCUT
 1293:             T(1) = T(1) + (MODRIJ - OPP_RCUT) * &
 1294:                 (C * OPP_K * (OPP_RCUT - 1) + S - &
 1295:                 S * 3.D0 * (OPP_RCUT - 1) / OPP_RCUT ) / OPP_RCUT**3
 1296:             ! Terms to make the second derivative go to zero at OPP_RCUT
 1297: !            T(1) = T(1) - 0.5D0 * (MODYIJ - OPP_RCUT)**2 * &
 1298: !                (S * OPP_K**2 * (OPP_RCUT - 1) / OPP_RCUT**3 - &
 1299: !                C * 2.D0 * OPP_K / OPP_RCUT**3 + &
 1300: !                C * 6.D0 * OPP_K * (OPP_RCUT - 1) / OPP_RCUT**4 + &
 1301: !                S * 6.D0 / OPP_RCUT**4 - &
 1302: !                S * 12.D0 * (OPP_RCUT - 1) / OPP_RCUT**5)
1331:         END IF1303:         END IF
1332:         IF (BTEST(OPP_PARAMS, 1)) THEN1304:         IF (BTEST(OPP_PARAMS, 1)) THEN
1333:             CALC_GRAD_PARAMS(2) = CALC_DVDPHI(MODRIJ) * CALC_XPLOR(MODRIJ)1305:             ! OPP_PHI derivative
 1306:             ! Unmodified potential term
 1307:             T(2) = SIN(OPP_K * (MODRIJ - 1) - OPP_PHI) / MODRIJ**3
 1308:             ! Term to make the potential go to zero at OPP_RCUT
 1309:             T(2) = T(2) - S / OPP_RCUT **3
 1310:             ! Terms to make the first derivative go to zero at OPP_RCUT
 1311:             T(2) = T(2) + (MODRIJ - OPP_RCUT) * &
 1312:                 ( - C * OPP_K + S * 3.D0 / OPP_RCUT ) / OPP_RCUT**3
 1313:             ! Terms to make the second derivative go to zero at OPP_RCUT
 1314: !            T(2) = T(2) - 0.5D0 * (MODYIJ - OPP_RCUT)**2 * &
 1315: !                ( - S * OPP_K**2 / OPP_RCUT**3 - &
 1316: !                C * 6.D0 * OPP_K / OPP_RCUT**4 + S * 12.D0 / OPP_RCUT**5)
1334:         END IF1317:         END IF
1335: 1318: 
 1319:         CALC_GRAD_PARAMS(:) = T(:)
 1320: 
1336:     END FUNCTION CALC_GRAD_PARAMS1321:     END FUNCTION CALC_GRAD_PARAMS
1337: 1322: 
1338: !-------------------------------------------------------------------------------1323: !-------------------------------------------------------------------------------
1339: !1324: !
1340: ! Calculates the second derivative contributions for the potential parameters1325: ! Calculates the second derivative contributions for the potential parameters
1341: ! for a cluster1326: ! for a cluster
1342: ! RIJ: vector between particles1327: ! RIJ: vector between particles
1343: ! IDI, IDJ: indices of the atom pair1328: ! IDI, IDJ: indices of the atom pair
1344: !1329: !
1345: !-------------------------------------------------------------------------------1330: !-------------------------------------------------------------------------------
1351: 1336: 
1352:         ! HESS: the Hessian matrix1337:         ! HESS: the Hessian matrix
1353:         USE MODHESS, ONLY: HESS1338:         USE MODHESS, ONLY: HESS
1354: 1339: 
1355:         IMPLICIT NONE1340:         IMPLICIT NONE
1356: 1341: 
1357:         DOUBLE PRECISION, INTENT(IN) :: RIJ(3)1342:         DOUBLE PRECISION, INTENT(IN) :: RIJ(3)
1358:         INTEGER, INTENT(IN) :: IDI, IDJ1343:         INTEGER, INTENT(IN) :: IDI, IDJ
1359: 1344: 
1360:         ! MODRIJ: particle-particle distance1345:         ! MODRIJ: particle-particle distance
1361:         ! TEMP_GRAD: temporary calculation store1346:         ! S, C, SC, CC: recurring sin and cos terms
1362:         DOUBLE PRECISION :: MODRIJ, TEMP_GRAD(3)1347:         ! TEMP, TEMP_GRAD: temporary calculation stores
 1348:         ! D2VDK2, D2VDPHI2, D2VDKDPHI: pure parameter second derivatives
 1349:         ! D2VDRDK, D2VDRDPHI; mixed second derivatives of pair potential
 1350:         DOUBLE PRECISION :: MODRIJ, S, C, SC, CC, TEMP, TEMP_GRAD(3)
 1351:         DOUBLE PRECISION :: D2VDK2, D2VDPHI2, D2VDKDPHI, D2VDRDK, D2VDRDPHI
 1352: 
 1353:         ! I: loop index
 1354:         INTEGER :: I
1363: 1355: 
1364:         ! Calculate some factors1356:         ! Calculate some factors
1365:         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))1357:         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))
 1358:         S = SIN(OPP_K * (MODRIJ - 1) - OPP_PHI)
 1359:         C = COS(OPP_K * (MODRIJ - 1) - OPP_PHI)
 1360:         SC = SIN(OPP_K * (OPP_RCUT - 1) - OPP_PHI)
 1361:         CC = COS(OPP_K * (OPP_RCUT - 1) - OPP_PHI)
1366: 1362: 
1367:         IF (BTEST(OPP_PARAMS, 0)) THEN1363:         IF (BTEST(OPP_PARAMS, 0)) THEN
1368:             ! OPP_K pure derivative1364:             ! OPP_K pure derivative
1369:             HESS(3*NATOMS-2,3*NATOMS-2) = HESS(3*NATOMS-2,3*NATOMS-2) + &1365:             ! Unmodified potential term
1370:                 CALC_D2VDK2(MODRIJ) * CALC_XPLOR(MODRIJ)1366:             D2VDK2 = - C * (MODRIJ - 1)**2 / MODRIJ**3
 1367:             HESS(3*NATOMS-2,3*NATOMS-2) = HESS(3*NATOMS-2,3*NATOMS-2) + D2VDK2
1371: 1368: 
1372:             ! OPP_K and position mixed derivative1369:             ! OPP_K and position mixed derivative
1373:             TEMP_GRAD(:) = (CALC_D2VDRDK(MODRIJ) * CALC_XPLOR(MODRIJ) + &1370:             ! Unmodified potential term
1374:                 CALC_DVDK(MODRIJ) * CALC_DXPLORDR(MODRIJ)) * RIJ(:) / MODRIJ1371:             D2VDRDK = (- C * OPP_K * (MODRIJ - 1) - S + &
 1372:                 S * 3.D0 * (MODRIJ - 1) / MODRIJ) / MODRIJ**3
 1373:             D2VDRDK = D2VDRDK / MODRIJ
 1374:             TEMP_GRAD(:) = D2VDRDK * RIJ(:)
1375:             HESS(3*IDI-2:3*IDI,3*NATOMS-2) = &1375:             HESS(3*IDI-2:3*IDI,3*NATOMS-2) = &
1376:                 HESS(3*IDI-2:3*IDI,3*NATOMS-2) + TEMP_GRAD(:)1376:                 HESS(3*IDI-2:3*IDI,3*NATOMS-2) + TEMP_GRAD(:)
1377:             HESS(3*NATOMS-2,3*IDI-2:3*IDI) = &1377:             HESS(3*NATOMS-2,3*IDI-2:3*IDI) = &
1378:                 HESS(3*NATOMS-2,3*IDI-2:3*IDI) + TEMP_GRAD(:)1378:                 HESS(3*NATOMS-2,3*IDI-2:3*IDI) + TEMP_GRAD(:)
1379:             HESS(3*IDJ-2:3*IDJ,3*NATOMS-2) = &1379:             HESS(3*IDJ-2:3*IDJ,3*NATOMS-2) = &
1380:                 HESS(3*IDJ-2:3*IDJ,3*NATOMS-2) - TEMP_GRAD(:)1380:                 HESS(3*IDJ-2:3*IDJ,3*NATOMS-2) - TEMP_GRAD(:)
1381:             HESS(3*NATOMS-2,3*IDJ-2:3*IDJ) = &1381:             HESS(3*NATOMS-2,3*IDJ-2:3*IDJ) = &
1382:                 HESS(3*NATOMS-2,3*IDJ-2:3*IDJ) - TEMP_GRAD(:)1382:                 HESS(3*NATOMS-2,3*IDJ-2:3*IDJ) - TEMP_GRAD(:)
1383:         END IF        1383:         END IF        
1384: 1384: 
1385:         IF (BTEST(OPP_PARAMS, 1)) THEN1385:         IF (BTEST(OPP_PARAMS, 1)) THEN
1386:             ! OPP_PHI pure derivative1386:             ! OPP_PHI pure derivative
1387:             HESS(3*NATOMS-1,3*NATOMS-1) = HESS(3*NATOMS-1,3*NATOMS-1) + &1387:             ! Unmodified potential term
1388:                 CALC_D2VDPHI2(MODRIJ) * CALC_XPLOR(MODRIJ)1388:             D2VDPHI2 = - C / MODRIJ**3
 1389:             HESS(3*NATOMS-1,3*NATOMS-1) = HESS(3*NATOMS-1,3*NATOMS-1) + D2VDPHI2
1389: 1390: 
1390:             ! OPP_PHI and position mixed derivative1391:             ! OPP_PHI and position mixed derivative
1391:             TEMP_GRAD(:) = (CALC_D2VDRDPHI(MODRIJ) * CALC_XPLOR(MODRIJ) + &1392:             ! Unmodified potential term
1392:                 CALC_DVDPHI(MODRIJ) * CALC_DXPLORDR(MODRIJ)) * RIJ(:) / MODRIJ1393:             D2VDRDPHI = (C * OPP_K - S * 3.D0 / MODRIJ) / MODRIJ**3
 1394:             D2VDRDPHI = D2VDRDPHI  / MODRIJ
 1395:             TEMP_GRAD(:) = D2VDRDPHI * RIJ(:)
1393:             HESS(3*IDI-2:3*IDI,3*NATOMS-1) = &1396:             HESS(3*IDI-2:3*IDI,3*NATOMS-1) = &
1394:                 HESS(3*IDI-2:3*IDI,3*NATOMS-1) + TEMP_GRAD(:)1397:                 HESS(3*IDI-2:3*IDI,3*NATOMS-1) + TEMP_GRAD(:)
1395:             HESS(3*NATOMS-1,3*IDI-2:3*IDI) = &1398:             HESS(3*NATOMS-1,3*IDI-2:3*IDI) = &
1396:                 HESS(3*NATOMS-1,3*IDI-2:3*IDI) + TEMP_GRAD(:)1399:                 HESS(3*NATOMS-1,3*IDI-2:3*IDI) + TEMP_GRAD(:)
1397:             HESS(3*IDJ-2:3*IDJ,3*NATOMS-1) = &1400:             HESS(3*IDJ-2:3*IDJ,3*NATOMS-1) = &
1398:                 HESS(3*IDJ-2:3*IDJ,3*NATOMS-1) - TEMP_GRAD(:)1401:                 HESS(3*IDJ-2:3*IDJ,3*NATOMS-1) - TEMP_GRAD(:)
1399:             HESS(3*NATOMS-1,3*IDJ-2:3*IDJ) = &1402:             HESS(3*NATOMS-1,3*IDJ-2:3*IDJ) = &
1400:                 HESS(3*NATOMS-1,3*IDJ-2:3*IDJ) - TEMP_GRAD(:)1403:                 HESS(3*NATOMS-1,3*IDJ-2:3*IDJ) - TEMP_GRAD(:)
1401:         END IF1404:         END IF
1402: 1405: 
1403:         IF (BTEST(OPP_PARAMS, 0) .AND. BTEST(OPP_PARAMS, 1)) THEN1406:         IF (BTEST(OPP_PARAMS, 0) .AND. BTEST(OPP_PARAMS, 1)) THEN
1404:             ! OPP_K and OPP_PHI derivative1407:             ! OPP_K and OPP_PHI derivative
1405:             HESS(3*NATOMS-2,3*NATOMS-1) = HESS(3*NATOMS-2,3*NATOMS-1) + &1408:             ! Unmodified potential term
1406:                 CALC_D2VDKDPHI(MODRIJ) * CALC_XPLOR(MODRIJ)1409:             D2VDKDPHI = C * (MODRIJ - 1) / MODRIJ**3
1407:             HESS(3*NATOMS-1,3*NATOMS-2) = HESS(3*NATOMS-1,3*NATOMS-2) + &1410:             HESS(3*NATOMS-2,3*NATOMS-1) = HESS(3*NATOMS-2,3*NATOMS-1)+D2VDKDPHI
1408:                 CALC_D2VDKDPHI(MODRIJ) * CALC_XPLOR(MODRIJ)1411:             HESS(3*NATOMS-1,3*NATOMS-2) = HESS(3*NATOMS-1,3*NATOMS-2)+D2VDKDPHI
1409:         END IF1412:         END IF
1410: 1413: 
1411:     END SUBROUTINE CALC_HESS_PARAMS1414:     END SUBROUTINE CALC_HESS_PARAMS
1412: 1415: 
1413: !-------------------------------------------------------------------------------1416: !-------------------------------------------------------------------------------
1414: !1417: !
1415: ! Calculates the second derivative contributions for the potential parameters in1418: ! Calculates the second derivative contributions for the potential parameters in
1416: ! a periodic system.1419: ! a periodic system.
1417: ! RIJ: vector between particles in fractional coordinates1420: ! RIJ: vector between particles in fractional coordinates
1418: ! BD: box derivative information1421: ! BD: box derivative information
1432:         USE TRICLINIC_MOD, ONLY: BOX_DERIV1435:         USE TRICLINIC_MOD, ONLY: BOX_DERIV
1433: 1436: 
1434:         IMPLICIT NONE1437:         IMPLICIT NONE
1435: 1438: 
1436:         DOUBLE PRECISION, INTENT(IN) :: RIJ(3)1439:         DOUBLE PRECISION, INTENT(IN) :: RIJ(3)
1437:         TYPE(BOX_DERIV), INTENT(IN) :: BD1440:         TYPE(BOX_DERIV), INTENT(IN) :: BD
1438:         INTEGER, INTENT(IN) :: IDI, IDJ1441:         INTEGER, INTENT(IN) :: IDI, IDJ
1439: 1442: 
1440:         ! YIJ: particle-particle vector in Cartesian space1443:         ! YIJ: particle-particle vector in Cartesian space
1441:         ! MODYIJ: particle-particle distance in Cartesian space1444:         ! MODYIJ: particle-particle distance in Cartesian space
 1445:         ! S, C, SC, CC: recurring sin and cos terms
1442:         ! TEMP, TEMP_GRAD: temporary calculation stores1446:         ! TEMP, TEMP_GRAD: temporary calculation stores
1443:         DOUBLE PRECISION :: YIJ(3), MODYIJ, TEMP, TEMP_GRAD(3)1447:         ! D2VDK2, D2VDPHI2, D2VDKDPHI: pure parameter second derivatives
 1448:         ! D2VDYDK, D2VDYDPHI; mixed second derivatives of pair potential
 1449:         DOUBLE PRECISION :: YIJ(3), MODYIJ, S, C, SC, CC, TEMP, TEMP_GRAD(3)
 1450:         DOUBLE PRECISION :: D2VDK2, D2VDPHI2, D2VDKDPHI, D2VDYDK, D2VDYDPHI
1444: 1451: 
1445:         ! I: loop index1452:         ! I: loop index
1446:         INTEGER :: I1453:         INTEGER :: I
1447: 1454: 
1448:         ! Calculate some factors1455:         ! Calculate some factors
1449:         YIJ(:) = MATMUL(BD%ORTHOG(:,:), RIJ(:))1456:         YIJ(:) = MATMUL(BD%ORTHOG(:,:), RIJ(:))
1450:         MODYIJ = SQRT(DOT_PRODUCT(YIJ(:), YIJ(:)))1457:         MODYIJ = SQRT(DOT_PRODUCT(YIJ(:), YIJ(:)))
 1458:         S = SIN(OPP_K * (MODYIJ - 1) - OPP_PHI)
 1459:         C = COS(OPP_K * (MODYIJ - 1) - OPP_PHI)
 1460:         SC = SIN(OPP_K * (OPP_RCUT - 1) - OPP_PHI)
 1461:         CC = COS(OPP_K * (OPP_RCUT - 1) - OPP_PHI)
1451: 1462: 
1452:         IF (BTEST(OPP_PARAMS, 0)) THEN1463:         IF (BTEST(OPP_PARAMS, 0)) THEN
1453:             ! OPP_K pure derivative1464:             ! OPP_K pure derivative
1454:             HESS(3*NATOMS-8,3*NATOMS-8) = HESS(3*NATOMS-8,3*NATOMS-8) + &1465:             ! Unmodified potential term
1455:                 CALC_D2VDK2(MODYIJ) * CALC_XPLOR(MODYIJ) * &1466:             D2VDK2 = - C * (MODYIJ - 1)**2 / MODYIJ**3
1456:                 MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ)1467:             ! Term to make the potential go to zero at OPP_RCUT
 1468:             D2VDK2 = D2VDK2 + CC * (OPP_RCUT - 1)**2 / OPP_RCUT **3
 1469:             ! Terms to make the first derivative go to zero at OPP_RCUT
 1470:             D2VDK2 = D2VDK2 + (MODYIJ - OPP_RCUT) * &
 1471:                 (- SC * OPP_K * (OPP_RCUT-1)**2 + CC * 2.D0 * (OPP_RCUT-1) - &
 1472:                 CC * 3.D0 * (OPP_RCUT - 1)**2 / OPP_RCUT ) / OPP_RCUT**3
 1473:             D2VDK2 = D2VDK2 * MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ)
 1474:             HESS(3*NATOMS-8,3*NATOMS-8) = HESS(3*NATOMS-8,3*NATOMS-8) + D2VDK2
1457: 1475: 
1458:             ! OPP_K and position mixed derivative1476:             ! OPP_K and position mixed derivative
1459:             TEMP_GRAD(:) = (CALC_D2VDRDK(MODYIJ) * CALC_XPLOR(MODYIJ) + &1477:             ! Unmodified potential term
1460:                 CALC_DVDK(MODYIJ) * CALC_DXPLORDR(MODYIJ)) * &1478:             D2VDYDK = (- C * OPP_K * (MODYIJ - 1) - S + &
1461:                 MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * &1479:                 S * 3.D0 * (MODYIJ - 1) / MODYIJ) / MODYIJ**3
1462:                 MATMUL(TRANSPOSE(BD%ORTHOG(:,:)), YIJ(:)) / MODYIJ1480:             ! Terms to make the first derivative go to zero at OPP_RCUT
 1481:             D2VDYDK = D2VDYDK + (CC * OPP_K * (OPP_RCUT - 1) + SC - &
 1482:                 SC * 3.D0 * (OPP_RCUT - 1) / OPP_RCUT) / OPP_RCUT**3
 1483:             D2VDYDK = D2VDYDK * MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) / MODYIJ
 1484:             TEMP_GRAD(:) = D2VDYDK * MATMUL(TRANSPOSE(BD%ORTHOG(:,:)), YIJ(:))
1463:             HESS(3*IDI-2:3*IDI,3*NATOMS-8) = &1485:             HESS(3*IDI-2:3*IDI,3*NATOMS-8) = &
1464:                 HESS(3*IDI-2:3*IDI,3*NATOMS-8) + TEMP_GRAD(:)1486:                 HESS(3*IDI-2:3*IDI,3*NATOMS-8) + TEMP_GRAD(:)
1465:             HESS(3*NATOMS-8,3*IDI-2:3*IDI) = &1487:             HESS(3*NATOMS-8,3*IDI-2:3*IDI) = &
1466:                 HESS(3*NATOMS-8,3*IDI-2:3*IDI) + TEMP_GRAD(:)1488:                 HESS(3*NATOMS-8,3*IDI-2:3*IDI) + TEMP_GRAD(:)
1467:             HESS(3*IDJ-2:3*IDJ,3*NATOMS-8) = &1489:             HESS(3*IDJ-2:3*IDJ,3*NATOMS-8) = &
1468:                 HESS(3*IDJ-2:3*IDJ,3*NATOMS-8) - TEMP_GRAD(:)1490:                 HESS(3*IDJ-2:3*IDJ,3*NATOMS-8) - TEMP_GRAD(:)
1469:             HESS(3*NATOMS-8,3*IDJ-2:3*IDJ) = &1491:             HESS(3*NATOMS-8,3*IDJ-2:3*IDJ) = &
1470:                 HESS(3*NATOMS-8,3*IDJ-2:3*IDJ) - TEMP_GRAD(:)1492:                 HESS(3*NATOMS-8,3*IDJ-2:3*IDJ) - TEMP_GRAD(:)
1471: 1493: 
1472:             ! OPP_K and unit cell parameter mixed derivative1494:             ! OPP_K and unit cell parameter mixed derivative
1473:             DO I = 1, 6 ! loop over unit cell parameters1495:             DO I = 1, 6 ! loop over unit cell parameters
1474:                 TEMP = (CALC_D2VDRDK(MODYIJ) * CALC_XPLOR(MODYIJ) + &1496:                 TEMP = D2VDYDK * &
1475:                     CALC_DVDK(MODYIJ) * CALC_DXPLORDR(MODYIJ)) * &1497:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,I), RIJ(:)))
1476:                     MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * & 
1477:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,I), RIJ(:))) / & 
1478:                     MODYIJ 
1479:                 HESS(3*NATOMS-6+I,3*NATOMS-8) = &1498:                 HESS(3*NATOMS-6+I,3*NATOMS-8) = &
1480:                     HESS(3*NATOMS-6+I,3*NATOMS-8) + TEMP1499:                     HESS(3*NATOMS-6+I,3*NATOMS-8) + TEMP
1481:                 HESS(3*NATOMS-8,3*NATOMS-6+I) = &1500:                 HESS(3*NATOMS-8,3*NATOMS-6+I) = &
1482:                     HESS(3*NATOMS-8,3*NATOMS-6+I) + TEMP1501:                     HESS(3*NATOMS-8,3*NATOMS-6+I) + TEMP
1483:             END DO ! end loop over unit cell parameters1502:             END DO ! end loop over unit cell parameters
1484:         END IF        1503:         END IF        
1485: 1504: 
1486:         IF (BTEST(OPP_PARAMS, 1)) THEN1505:         IF (BTEST(OPP_PARAMS, 1)) THEN
1487:             ! OPP_PHI pure derivative1506:             ! OPP_PHI pure derivative
1488:             HESS(3*NATOMS-7,3*NATOMS-7) = HESS(3*NATOMS-7,3*NATOMS-7) + &1507:             ! Unmodified potential term
1489:                 CALC_D2VDPHI2(MODYIJ) * CALC_XPLOR(MODYIJ) * &1508:             D2VDPHI2 = - C / MODYIJ**3
1490:                 MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ)1509:             ! Term to make the potential go to zero at OPP_RCUT
 1510:             D2VDPHI2 = D2VDPHI2 + CC / OPP_RCUT **3
 1511:             ! Terms to make the first derivative go to zero at OPP_RCUT
 1512:             D2VDPHI2 = D2VDPHI2 + (MODYIJ - OPP_RCUT) * &
 1513:                 (- SC * OPP_K - CC * 3.D0 / OPP_RCUT ) / OPP_RCUT**3
 1514:             D2VDPHI2 = D2VDPHI2 * MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ)
 1515:             HESS(3*NATOMS-7,3*NATOMS-7) = HESS(3*NATOMS-7,3*NATOMS-7) + D2VDPHI2
1491: 1516: 
1492:             ! OPP_PHI and position mixed derivative1517:             ! OPP_PHI and position mixed derivative
1493:             TEMP_GRAD(:) = (CALC_D2VDRDPHI(MODYIJ) * CALC_XPLOR(MODYIJ) + &1518:             ! Unmodified potential term
1494:                 CALC_DVDPHI(MODYIJ) * CALC_DXPLORDR(MODYIJ)) * &1519:             D2VDYDPHI = (C * OPP_K - S * 3.D0 / MODYIJ) / MODYIJ**3
1495:                 MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * &1520:             ! Terms to make the first derivative go to zero at OPP_RCUT
1496:                 MATMUL(TRANSPOSE(BD%ORTHOG(:,:)), YIJ(:)) / MODYIJ1521:             D2VDYDPHI = D2VDYDPHI + ( - CC * OPP_K + &
 1522:                 SC * 3.D0 / OPP_RCUT) / OPP_RCUT**3
 1523:             D2VDYDPHI = D2VDYDPHI * MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) / MODYIJ
 1524:             TEMP_GRAD(:) = D2VDYDPHI * MATMUL(TRANSPOSE(BD%ORTHOG(:,:)), YIJ(:))
1497:             HESS(3*IDI-2:3*IDI,3*NATOMS-7) = &1525:             HESS(3*IDI-2:3*IDI,3*NATOMS-7) = &
1498:                 HESS(3*IDI-2:3*IDI,3*NATOMS-7) + TEMP_GRAD(:)1526:                 HESS(3*IDI-2:3*IDI,3*NATOMS-7) + TEMP_GRAD(:)
1499:             HESS(3*NATOMS-7,3*IDI-2:3*IDI) = &1527:             HESS(3*NATOMS-7,3*IDI-2:3*IDI) = &
1500:                 HESS(3*NATOMS-7,3*IDI-2:3*IDI) + TEMP_GRAD(:)1528:                 HESS(3*NATOMS-7,3*IDI-2:3*IDI) + TEMP_GRAD(:)
1501:             HESS(3*IDJ-2:3*IDJ,3*NATOMS-7) = &1529:             HESS(3*IDJ-2:3*IDJ,3*NATOMS-7) = &
1502:                 HESS(3*IDJ-2:3*IDJ,3*NATOMS-7) - TEMP_GRAD(:)1530:                 HESS(3*IDJ-2:3*IDJ,3*NATOMS-7) - TEMP_GRAD(:)
1503:             HESS(3*NATOMS-7,3*IDJ-2:3*IDJ) = &1531:             HESS(3*NATOMS-7,3*IDJ-2:3*IDJ) = &
1504:                 HESS(3*NATOMS-7,3*IDJ-2:3*IDJ) - TEMP_GRAD(:)1532:                 HESS(3*NATOMS-7,3*IDJ-2:3*IDJ) - TEMP_GRAD(:)
1505: 1533: 
1506:             ! OPP_K and unit cell parameter mixed derivative1534:             ! OPP_K and unit cell parameter mixed derivative
1507:             DO I = 1, 6 ! loop over unit cell parameters1535:             DO I = 1, 6 ! loop over unit cell parameters
1508:                 TEMP = (CALC_D2VDRDPHI(MODYIJ) * CALC_XPLOR(MODYIJ) + &1536:                 TEMP = D2VDYDPHI * &
1509:                     CALC_DVDPHI(MODYIJ) * CALC_DXPLORDR(MODYIJ)) * &1537:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,I), RIJ(:)))
1510:                     MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * & 
1511:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,I), RIJ(:))) / & 
1512:                     MODYIJ 
1513:                 HESS(3*NATOMS-6+I,3*NATOMS-7) = &1538:                 HESS(3*NATOMS-6+I,3*NATOMS-7) = &
1514:                     HESS(3*NATOMS-6+I,3*NATOMS-7) + TEMP1539:                     HESS(3*NATOMS-6+I,3*NATOMS-7) + TEMP
1515:                 HESS(3*NATOMS-7,3*NATOMS-6+I) = &1540:                 HESS(3*NATOMS-7,3*NATOMS-6+I) = &
1516:                     HESS(3*NATOMS-7,3*NATOMS-6+I) + TEMP1541:                     HESS(3*NATOMS-7,3*NATOMS-6+I) + TEMP
1517:             END DO ! end loop over unit cell parameters1542:             END DO ! end loop over unit cell parameters
1518:         END IF1543:         END IF
1519: 1544: 
1520:         IF (BTEST(OPP_PARAMS, 0) .AND. BTEST(OPP_PARAMS, 1)) THEN1545:         IF (BTEST(OPP_PARAMS, 0) .AND. BTEST(OPP_PARAMS, 1)) THEN
1521:             ! OPP_K and OPP_PHI derivative1546:             ! OPP_K and OPP_PHI derivative
1522:             TEMP = CALC_D2VDKDPHI(MODYIJ) * CALC_XPLOR(MODYIJ) * &1547:             ! Unmodified potential term
1523:                 MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ)1548:             D2VDKDPHI = C * (MODYIJ - 1) / MODYIJ**3
1524:             HESS(3*NATOMS-8,3*NATOMS-7) = HESS(3*NATOMS-8,3*NATOMS-7) + TEMP1549:             ! Term to make the potential go to zero at OPP_RCUT
1525:             HESS(3*NATOMS-7,3*NATOMS-8) = HESS(3*NATOMS-7,3*NATOMS-8) + TEMP1550:             D2VDKDPHI = D2VDKDPHI - CC * (OPP_RCUT - 1) / OPP_RCUT**3
 1551:             ! Terms to make the first derivative go to zero at OPP_RCUT
 1552:             D2VDKDPHI = D2VDKDPHI + (MODYIJ - OPP_RCUT) * &
 1553:                 (SC * OPP_K * (OPP_RCUT - 1) -  CC + &
 1554:                 CC * 3.D0 * (OPP_RCUT - 1) / OPP_RCUT) / OPP_RCUT**3
 1555:             D2VDKDPHI = D2VDKDPHI * MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ)
 1556:             HESS(3*NATOMS-8,3*NATOMS-7) = HESS(3*NATOMS-8,3*NATOMS-7)+D2VDKDPHI
 1557:             HESS(3*NATOMS-7,3*NATOMS-8) = HESS(3*NATOMS-7,3*NATOMS-8)+D2VDKDPHI
1526:         END IF1558:         END IF
1527: 1559: 
1528:     END SUBROUTINE CALC_HESS_PARAMS_PER1560:     END SUBROUTINE CALC_HESS_PARAMS_PER
1529: 1561: 
1530: !-------------------------------------------------------------------------------1562: !-------------------------------------------------------------------------------
1531: !1563: !
1532: ! Derivative of potential wrt OPP_K without XPLOR 
1533: ! MODRIJ: distance between particles 
1534: ! 
1535: !------------------------------------------------------------------------------- 
1536:     PURE DOUBLE PRECISION FUNCTION CALC_DVDK (MODRIJ) 
1537:  
1538:         IMPLICIT NONE 
1539:  
1540:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
1541:  
1542:         CALC_DVDK = - SIN(OPP_K * (MODRIJ-1) - OPP_PHI) * (MODRIJ-1) / MODRIJ**3 
1543:  
1544:     END FUNCTION CALC_DVDK 
1545:  
1546: !------------------------------------------------------------------------------- 
1547: ! 
1548: ! Derivative of potential wrt OPP_PHI without XPLOR 
1549: ! MODRIJ: distance between particles 
1550: ! 
1551: !------------------------------------------------------------------------------- 
1552:     PURE DOUBLE PRECISION FUNCTION CALC_DVDPHI (MODRIJ) 
1553:  
1554:         IMPLICIT NONE 
1555:  
1556:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
1557:  
1558:         CALC_DVDPHI = SIN(OPP_K * (MODRIJ - 1) - OPP_PHI) / MODRIJ**3 
1559:  
1560:     END FUNCTION CALC_DVDPHI 
1561:  
1562: !------------------------------------------------------------------------------- 
1563: ! 
1564: ! Second derivative of potential wrt OPP_K twice without XPLOR 
1565: ! MODRIJ: distance between particles 
1566: ! 
1567: !------------------------------------------------------------------------------- 
1568:     PURE DOUBLE PRECISION FUNCTION CALC_D2VDK2 (MODRIJ) 
1569:  
1570:         IMPLICIT NONE 
1571:  
1572:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
1573:  
1574:         CALC_D2VDK2 = & 
1575:             - COS(OPP_K * (MODRIJ - 1) - OPP_PHI) * (MODRIJ - 1)**2 / MODRIJ**3 
1576:  
1577:     END FUNCTION CALC_D2VDK2 
1578:  
1579: !------------------------------------------------------------------------------- 
1580: ! 
1581: ! Second derivative of potential wrt OPP_PHI twice without XPLOR 
1582: ! MODRIJ: distance between particles 
1583: ! 
1584: !------------------------------------------------------------------------------- 
1585:     PURE DOUBLE PRECISION FUNCTION CALC_D2VDPHI2 (MODRIJ) 
1586:  
1587:         IMPLICIT NONE 
1588:  
1589:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
1590:  
1591:         CALC_D2VDPHI2 = - COS(OPP_K * (MODRIJ - 1) - OPP_PHI) / MODRIJ**3 
1592:  
1593:     END FUNCTION CALC_D2VDPHI2 
1594:  
1595: !------------------------------------------------------------------------------- 
1596: ! 
1597: ! Second derivative of potential wrt OPP_K and OPP_PHI without XPLOR 
1598: ! MODRIJ: distance between particles 
1599: ! 
1600: !------------------------------------------------------------------------------- 
1601:     PURE DOUBLE PRECISION FUNCTION CALC_D2VDKDPHI (MODRIJ) 
1602:  
1603:         IMPLICIT NONE 
1604:  
1605:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
1606:  
1607:         CALC_D2VDKDPHI = & 
1608:             COS(OPP_K * (MODRIJ - 1) - OPP_PHI) * (MODRIJ - 1) / MODRIJ**3 
1609:  
1610:     END FUNCTION CALC_D2VDKDPHI 
1611:  
1612: !------------------------------------------------------------------------------- 
1613: ! 
1614: ! Second derivative of potential wrt position and OPP_K without XPLOR 
1615: ! MODRIJ: distance between particles 
1616: ! 
1617: !------------------------------------------------------------------------------- 
1618:     PURE DOUBLE PRECISION FUNCTION CALC_D2VDRDK (MODRIJ) 
1619:  
1620:         IMPLICIT NONE 
1621:  
1622:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
1623:  
1624:         CALC_D2VDRDK = & 
1625:             (- COS(OPP_K * (MODRIJ - 1) - OPP_PHI) * OPP_K * (MODRIJ - 1) - & 
1626:             SIN(OPP_K * (MODRIJ - 1) - OPP_PHI) + & 
1627:             SIN(OPP_K * (MODRIJ - 1) - OPP_PHI) * 3.D0 * (MODRIJ - 1) / & 
1628:             MODRIJ) / MODRIJ**3 
1629:  
1630:     END FUNCTION CALC_D2VDRDK 
1631:  
1632: !------------------------------------------------------------------------------- 
1633: ! 
1634: ! Second derivative of potential wrt position and OPP_PHI without XPLOR 
1635: ! MODRIJ: distance between particles 
1636: ! 
1637: !------------------------------------------------------------------------------- 
1638:     PURE DOUBLE PRECISION FUNCTION CALC_D2VDRDPHI (MODRIJ) 
1639:  
1640:         IMPLICIT NONE 
1641:  
1642:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
1643:  
1644:         CALC_D2VDRDPHI = (COS(OPP_K * (MODRIJ - 1) - OPP_PHI) * OPP_K - & 
1645:             SIN(OPP_K * (MODRIJ - 1) - OPP_PHI) * 3.D0 / MODRIJ) / MODRIJ**3 
1646:  
1647:     END FUNCTION CALC_D2VDRDPHI 
1648:  
1649: !------------------------------------------------------------------------------- 
1650: ! 
1651: ! Repel the unit cell volume from approaching 0, which repels the unit cell from1564: ! Repel the unit cell volume from approaching 0, which repels the unit cell from
1652: ! adopting physically impossible angle combinations. Uses a WCA potential1565: ! adopting physically impossible angle combinations. Uses a WCA potential
1653: ! ENERGY: energy of the system1566: ! ENERGY: energy of the system
1654: ! GRAD_ANGLES: gradients for the box angles1567: ! GRAD_ANGLES: gradients for the box angles
1655: ! BD: information on the box and its derivatives1568: ! BD: information on the box and its derivatives
1656: ! GTEST: whether to calculate gradients1569: ! GTEST: whether to calculate gradients
1657: ! STEST: whether to calculate second derivatives1570: ! STEST: whether to calculate second derivatives
1658: !1571: !
1659: !-------------------------------------------------------------------------------1572: !-------------------------------------------------------------------------------
1660: 1573: 
1694:                 GRAD_ANGLES(:) = GRAD_ANGLES(:) +  DEDV * BD%V_DERIV(:)1607:                 GRAD_ANGLES(:) = GRAD_ANGLES(:) +  DEDV * BD%V_DERIV(:)
1695:             END IF1608:             END IF
1696: 1609: 
1697:             IF (STEST) THEN1610:             IF (STEST) THEN
1698:                 ! Add the Hessian contributions, if necessary1611:                 ! Add the Hessian contributions, if necessary
1699:                 DEDV = 24.D0 * V_EPS / V_SIGMA * &1612:                 DEDV = 24.D0 * V_EPS / V_SIGMA * &
1700:                     ((V_SIGMA / BD%V)**(7) - 2.D0 * (V_SIGMA / BD%V)**(13))1613:                     ((V_SIGMA / BD%V)**(7) - 2.D0 * (V_SIGMA / BD%V)**(13))
1701:                 D2EDV2 = 24.D0 * V_EPS / V_SIGMA**2 * &1614:                 D2EDV2 = 24.D0 * V_EPS / V_SIGMA**2 * &
1702:                     (26.D0 * (V_SIGMA/BD%V)**(14) - 7.D0 * (V_SIGMA/BD%V)**(8))1615:                     (26.D0 * (V_SIGMA/BD%V)**(14) - 7.D0 * (V_SIGMA/BD%V)**(8))
1703: 1616: 
1704:                 DO I = 1, 3 ! loop over angles1617:                 DO I = 1, 3 ! outer loop over angles
1705:                     HESS(3*NATOMS-5+I,3*NATOMS-5:3*NATOMS-3) = &1618:                     HESS(3*NATOMS-5+I,3*NATOMS-5:3*NATOMS-3) = &
1706:                         HESS(3*NATOMS-5+I,3*NATOMS-5:3*NATOMS-3) + &1619:                         HESS(3*NATOMS-5+I,3*NATOMS-5:3*NATOMS-3) + &
1707:                         D2EDV2 * BD%V_DERIV(I) * BD%V_DERIV(:) + &1620:                         D2EDV2 * BD%V_DERIV(I) * BD%V_DERIV(:) + &
1708:                         DEDV * BD%V_DERIV2(I, :)1621:                         DEDV * BD%V_DERIV2(I, :)
1709:                 END DO ! end loop over angles1622:                 END DO ! end outer loop over angles
1710:             END IF1623:             END IF
1711:         END IF1624:         END IF
1712: 1625: 
1713:     END SUBROUTINE CONSTRAIN_VOLUME1626:     END SUBROUTINE CONSTRAIN_VOLUME
1714: 1627: 
1715: !-------------------------------------------------------------------------------1628: !-------------------------------------------------------------------------------
1716: !1629: !
1717: ! Impose limits on the allowed range of the potential parameters OPP_K by adding1630: ! Impose limits on the allowed range of the potential parameters OPP_K by adding
1718: ! a harmonic repulsion outside the allowed range.1631: ! a harmonic repulsion outside the allowed range.
1719: ! ENERGY: energy of the system1632: ! ENERGY: energy of the system
1799:         ! QMINP: coordinates of the saved minima1712:         ! QMINP: coordinates of the saved minima
1800:         ! FF: step at which the saved minima was first found1713:         ! FF: step at which the saved minima was first found
1801:         USE QMODULE, ONLY: QMIN, QMINP, FF1714:         USE QMODULE, ONLY: QMIN, QMINP, FF
1802: 1715: 
1803:         ! BOX_DERIV: type for storing box derivative information1716:         ! BOX_DERIV: type for storing box derivative information
1804:         ! CALC_BOX_DERIV: calculates box derivative information1717:         ! CALC_BOX_DERIV: calculates box derivative information
1805:         USE TRICLINIC_MOD, ONLY: BOX_DERIV, CALC_BOX_DERIV1718:         USE TRICLINIC_MOD, ONLY: BOX_DERIV, CALC_BOX_DERIV
1806: 1719: 
1807:         IMPLICIT NONE1720:         IMPLICIT NONE
1808: 1721: 
1809:         ! GETUNIT: function for finding a free file unit1722:         ! GETUNIT: function for fiding a free file unit
1810:         ! FUNIT: output file unit1723:         ! FUNIT: output file unit
1811:         ! NMOL: actual number of particles1724:         ! NMOL: actual number of particles
1812:         ! I, J: loop iterators1725:         ! I, J: loop iterators
1813:         INTEGER :: GETUNIT, FUNIT, NMOL, I, J1726:         INTEGER :: GETUNIT, FUNIT, NMOL, I, J
1814: 1727: 
1815:         ! BOX_SIZE: convenience array for the box sizes1728:         ! BOX_SIZE: convenience array for the box sizes
1816:         ! BOX_ANGLES: convenience array for the box angles1729:         ! BOX_ANGLES: convenience array for the box angles
1817:         ! COORD: coordinate in orthogonal space1730:         ! COORD: coordinate in orthogonal space
1818:         DOUBLE PRECISION :: BOX_SIZE(3), BOX_ANGLES(3), COORD(3)1731:         DOUBLE PRECISION :: BOX_SIZE(3), BOX_ANGLES(3), COORD(3)
1819: 1732: 


r33470/perc.f90 2017-11-13 18:30:13.669646537 +0000 r33469/perc.f90 2017-11-13 18:30:14.785661300 +0000
 10: !   GMIN is distributed in the hope that it will be useful, 10: !   GMIN is distributed in the hope that it will be useful,
 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: SUBROUTINE PERC(P,NATOMS,PERCCUT,PERCT,DEBUG,MYUNIT,RIGID) 19: SUBROUTINE PERC(P,NATOMS,PERCCUT,PERCT,DEBUG,MYUNIT,RIGID)
 20:   USE COMMONS, ONLY : MODULARCURRENTN, MODULART, OPPT 20:   USE COMMONS, ONLY : MODULARCURRENTN,MODULART
 21:   USE OPP_MOD, ONLY: OPP_MODE 
 22:   IMPLICIT NONE 21:   IMPLICIT NONE
 23:  22: 
 24:   INTEGER, INTENT(IN)  :: NATOMS, MYUNIT 23:   INTEGER, INTENT(IN)  :: NATOMS, MYUNIT
 25:   LOGICAL, INTENT(IN)  :: DEBUG, RIGID 24:   LOGICAL, INTENT(IN)  :: DEBUG, RIGID
 26:   LOGICAL, INTENT(OUT) :: PERCT 25:   LOGICAL, INTENT(OUT) :: PERCT
 27:   DOUBLE PRECISION, INTENT(IN)  :: P(3*NATOMS), PERCCUT 26:   DOUBLE PRECISION, INTENT(IN)  :: P(3*NATOMS), PERCCUT
 28:   INTEGER NSITES, J1, J2, NCYCLE, DMIN1, DMAX1, NUNCON1 27:   INTEGER NSITES, J1, J2, NCYCLE, DMIN1, DMAX1, NUNCON1
 29:   INTEGER, ALLOCATABLE :: NDIST1(:) 28:   INTEGER, ALLOCATABLE :: NDIST1(:)
 30:   DOUBLE PRECISION DUMMY 29:   DOUBLE PRECISION DUMMY
 31:   LOGICAL CHANGED 30:   LOGICAL CHANGED
 32:   LOGICAL, ALLOCATABLE :: CON(:,:) 31:   LOGICAL, ALLOCATABLE :: CON(:,:)
 33:  32: 
 34:   IF(MODULART) THEN 33:   IF(MODULART) THEN
 35:     NSITES = MODULARCURRENTN  34:     NSITES = MODULARCURRENTN 
 36:   ELSE IF (OPPT) THEN 
 37: ! jwrm2> Need to remove alchemical and unit cell degrees of freedom. 
 38:     SELECT CASE (OPP_MODE) 
 39:     CASE(1) 
 40:       NSITES = NATOMS - 2 
 41:     CASE(2) 
 42:       NSITES = NATOMS - 3 
 43:     CASE(3) 
 44:       NSITES = NATOMS - 1 
 45:     END SELECT 
 46:   ELSE 35:   ELSE
 47:     NSITES = NATOMS 36:     NSITES = NATOMS
 48:   END IF 37:   END IF
 49:   IF (RIGID) THEN 38:   IF (RIGID) THEN
 50:     NSITES = NSITES/2 39:     NSITES = NSITES/2
 51:   ENDIF 40:   ENDIF
 52: ! WRITE(MYUNIT,'(A,2L5,3I6)') 'MODULART,RIGID,NSITES,NATOMS,MODULARCURRENTN=',MODULART,RIGID,NSITES,NATOMS,MODULARCURRENTN 41: ! WRITE(MYUNIT,'(A,2L5,3I6)') 'MODULART,RIGID,NSITES,NATOMS,MODULARCURRENTN=',MODULART,RIGID,NSITES,NATOMS,MODULARCURRENTN
 53:   IF(ALLOCATED(NDIST1)) DEALLOCATE(NDIST1) 42:   IF(ALLOCATED(NDIST1)) DEALLOCATE(NDIST1)
 54:   ALLOCATE(NDIST1(NSITES), CON(NSITES,NSITES)) 43:   ALLOCATE(NDIST1(NSITES), CON(NSITES,NSITES))
 55:  44: 
120: 109: 
121:   USE COMMONS, ONLY : atomingroup, modulart, modularcurrentn110:   USE COMMONS, ONLY : atomingroup, modulart, modularcurrentn
122:   IMPLICIT NONE111:   IMPLICIT NONE
123: 112: 
124:   INTEGER, INTENT(IN)  :: NATOMS, MYUNIT113:   INTEGER, INTENT(IN)  :: NATOMS, MYUNIT
125:   LOGICAL, INTENT(IN)  :: DEBUG, RIGID114:   LOGICAL, INTENT(IN)  :: DEBUG, RIGID
126:   DOUBLE PRECISION, INTENT(IN)  :: PERCCUT115:   DOUBLE PRECISION, INTENT(IN)  :: PERCCUT
127:   INTEGER NSITES, J1, J2, NCYCLE, DMIN1, DMAX1, NUNCON1, groupindex, nsitesskip,nunconnsave,k1,k2116:   INTEGER NSITES, J1, J2, NCYCLE, DMIN1, DMAX1, NUNCON1, groupindex, nsitesskip,nunconnsave,k1,k2
128:   INTEGER, ALLOCATABLE :: NDIST1(:)117:   INTEGER, ALLOCATABLE :: NDIST1(:)
129:   DOUBLE PRECISION DUMMY, COORDSNEW(3*NATOMS), P(3*NATOMS) !P is now intent inout118:   DOUBLE PRECISION DUMMY, COORDSNEW(3*NATOMS), P(3*NATOMS) !P is now intent inout
130:   LOGICAL CHANGED119:   LOGICAL CHANGED, PERCT
131:   LOGICAL, ALLOCATABLE :: CON(:,:)120:   LOGICAL, ALLOCATABLE :: CON(:,:)
132:  121:  
133: 122: 
134:   IF(MODULART) THEN123:   IF(MODULART) THEN
135:     NSITES = MODULARCURRENTN 124:     NSITES = MODULARCURRENTN 
136:   ELSE125:   ELSE
137:     NSITES = NATOMS126:     NSITES = NATOMS
138:   END IF127:   END IF
139:   IF (RIGID) THEN128:   IF (RIGID) THEN
140:     NSITES = NSITES/2129:     NSITES = NSITES/2
156: ! to check percolation only for these. Array P is being reordered so that each 145: ! to check percolation only for these. Array P is being reordered so that each 
157: ! detected disconnected group is populating it sequentially. We take always AMERI...khm..CONNECTED particles first.146: ! detected disconnected group is populating it sequentially. We take always AMERI...khm..CONNECTED particles first.
158:  groupindex = groupindex + 1147:  groupindex = groupindex + 1
159:  NSITESSKIP=NSITES-NUNCONNSAVE !this is 0 for the first loop148:  NSITESSKIP=NSITES-NUNCONNSAVE !this is 0 for the first loop
160:   IF(DEBUG) WRITE(MYUNIT,*) 'percgroup> Building group ', groupindex, nunconnsave, nsitesskip149:   IF(DEBUG) WRITE(MYUNIT,*) 'percgroup> Building group ', groupindex, nunconnsave, nsitesskip
161: ! NSITES=NUNCONNSAVE150: ! NSITES=NUNCONNSAVE
162: 151: 
163:   CON(1:NSITES,1:NSITES)=.FALSE.152:   CON(1:NSITES,1:NSITES)=.FALSE.
164:   DO J1=NSITESSKIP+1,NSITES153:   DO J1=NSITESSKIP+1,NSITES
165:     DO J2=J1+1,NSITES154:     DO J2=J1+1,NSITES
166:       DUMMY=(COORDSNEW(3*(J2-1)+1)-COORDSNEW(3*(J1-1)+1))**2+&155:       DUMMY=(COORDSNEW(3*(J2-1)+1)-COORDSNEW(3*(J1-1)+1))**2+(COORDSNEW(3*(J2-1)+2)-COORDSNEW(3*(J1-1)+2))**2+(COORDSNEW(3*(J2-1)+3)-COORDSNEW(3*(J1-1)+3))**2
167:             (COORDSNEW(3*(J2-1)+2)-COORDSNEW(3*(J1-1)+2))**2+& 
168:             (COORDSNEW(3*(J2-1)+3)-COORDSNEW(3*(J1-1)+3))**2 
169:       IF (DUMMY.LT.PERCCUT) THEN156:       IF (DUMMY.LT.PERCCUT) THEN
170:         CON(J2,J1)=.TRUE.157:         CON(J2,J1)=.TRUE.
171:         CON(J1,J2)=.TRUE.158:         CON(J1,J2)=.TRUE.
172:       ENDIF159:       ENDIF
173:     ENDDO160:     ENDDO
174:   ENDDO161:   ENDDO
175: 162: 
176: ! 163: ! 
177: ! Check that we have a percolating constraint network.164: ! Check that we have a percolating constraint network.
178: !165: !
199:       ENDIF186:       ENDIF
200:     ENDDO187:     ENDDO
201:     IF ((NDIST1(J1).GT.DMAX1).AND.(NDIST1(J1).NE.1000000)) DMAX1=NDIST1(J1)188:     IF ((NDIST1(J1).GT.DMAX1).AND.(NDIST1(J1).NE.1000000)) DMAX1=NDIST1(J1)
202:     IF (NDIST1(J1).LT.DMIN1) DMIN1=NDIST1(J1)189:     IF (NDIST1(J1).LT.DMIN1) DMIN1=NDIST1(J1)
203:     IF (NDIST1(J1).EQ.1000000) NUNCON1=NUNCON1+1190:     IF (NDIST1(J1).EQ.1000000) NUNCON1=NUNCON1+1
204:   ENDDO191:   ENDDO
205: !  PRINT *,'DMIN1,DMAX1,NUNCON1,NCYCLE,CHANGED=',DMIN1,DMAX1,NUNCON1,NCYCLE,CHANGED192: !  PRINT *,'DMIN1,DMAX1,NUNCON1,NCYCLE,CHANGED=',DMIN1,DMAX1,NUNCON1,NCYCLE,CHANGED
206:   IF (CHANGED) GOTO 55193:   IF (CHANGED) GOTO 55
207:   IF (DEBUG) WRITE(MYUNIT,'(3(A,I8))') ' percgroup> steps to first atom in the group converged in ',NCYCLE-1, &194:   IF (DEBUG) WRITE(MYUNIT,'(3(A,I8))') ' percgroup> steps to first atom in the group converged in ',NCYCLE-1, &
208:      &                    ' cycles; maximum=',DMAX1,' disconnected=',NUNCON1195:      &                    ' cycles; maximum=',DMAX1,' disconnected=',NUNCON1
 196:   PERCT=.TRUE.
209:   NUNCONNSAVE=NUNCON1197:   NUNCONNSAVE=NUNCON1
210:   IF (NUNCON1.GT.0) THEN198:   IF (NUNCON1.GT.0) THEN
211: ! reorder the coordinates199: ! reorder the coordinates
212:      K1=NSITESSKIP200:      K1=NSITESSKIP
213:      K2=0201:      K2=0
214: !        WRITE(MYUNIT,*) 'before reordering', k1, k2, nsitesskip,nuncon1202: !        WRITE(MYUNIT,*) 'before reordering', k1, k2, nsitesskip,nuncon1
215: !        WRITE(MYUNIT,*) NDIST1(:)203: !        WRITE(MYUNIT,*) NDIST1(:)
216: !        WRITE(MYUNIT,*) 'original coordinates: '204: !        WRITE(MYUNIT,*) 'original coordinates: '
217: !        WRITE(MYUNIT,*) COORDSNEW(:)205: !        WRITE(MYUNIT,*) COORDSNEW(:)
218:      DO J1=NSITESSKIP+1,NSITES206:      DO J1=NSITESSKIP+1,NSITES
219: !        WRITE(MYUNIT,*) 'sf344 debug> J1=', J1, NDIST1(J1)207: !        WRITE(MYUNIT,*) 'sf344 debug> J1=', J1, NDIST1(J1)
220:         IF(NDIST1(J1)<1000000) THEN208:         IF(NDIST1(J1)<1000000) THEN
221:           K1=K1+1209:           K1=K1+1
222:           atomingroup(k1)=groupindex210:           atomingroup(k1)=groupindex
223:           P(3*K1-2:3*K1)=COORDSNEW(3*J1-2:3*J1) !copy over coordinates for connected group211:           P(3*K1-2:3*K1)=COORDSNEW(3*J1-2:3*J1) !copy over coordinates for connected group
224:           IF(RIGID) THEN212:           IF(RIGID) P(3*K1-2+3*NATOMS/2:3*K1+3*NATOMS/2)=COORDSNEW(3*J1-2+3*NATOMS/2:3*J1+3*NATOMS/2) !copy over angle-axis coordinates as well
225:             !copy over angle-axis coordinates as well 
226:             P(3*K1-2+3*NATOMS/2:3*K1+3*NATOMS/2)=COORDSNEW(3*J1-2+3*NATOMS/2:3*J1+3*NATOMS/2) 
227:           END IF 
228:         ELSE !disconnected particle, coordinates should go to the end of the array213:         ELSE !disconnected particle, coordinates should go to the end of the array
229:           K2=K2+1214:           K2=K2+1
230: !          WRITE(MYUNIT,*) 'sf344 debug> ', 3*NSITES-3*NUNCON1+3*K2-2215: !          WRITE(MYUNIT,*) 'sf344 debug> ', 3*NSITES-3*NUNCON1+3*K2-2
231: !          WRITE(MYUNIT,*) COORDSNEW(3*J1-2:3*J1)216: !          WRITE(MYUNIT,*) COORDSNEW(3*J1-2:3*J1)
232:           P(3*NSITES-3*NUNCON1+3*K2-2:3*NSITES-3*NUNCON1+3*K2)=COORDSNEW(3*J1-2:3*J1)217:           P(3*NSITES-3*NUNCON1+3*K2-2:3*NSITES-3*NUNCON1+3*K2)=COORDSNEW(3*J1-2:3*J1)
233:           IF(RIGID) THEN218:          IF(RIGID) P(3*NSITES-3*NUNCON1+3*K2+3*NATOMS/2-2:3*NSITES-3*NUNCON1+3*K2+3*NATOMS/2)=COORDSNEW(3*J1-2+3*NATOMS/2:3*J1+3*NATOMS/2)
234:             P(3*NSITES-3*NUNCON1+3*K2+3*NATOMS/2-2:3*NSITES-3*NUNCON1+3*K2+3*NATOMS/2)=COORDSNEW(3*J1-2+3*NATOMS/2:3*J1+3*NATOMS/2) 
235:           END IF 
236:         END IF219:         END IF
237:      END DO220:      END DO
238: !        WRITE(MYUNIT,*) 'after reordering', k1, k2, nsitesskip221: !        WRITE(MYUNIT,*) 'after reordering', k1, k2, nsitesskip
239: !        WRITE(MYUNIT,*) 'reordered coordinates: '222: !        WRITE(MYUNIT,*) 'reordered coordinates: '
240: !        WRITE(MYUNIT,*) P(:)223: !        WRITE(MYUNIT,*) P(:)
241: ! save the reordered coordinate array224: ! save the reordered coordinate array
242:      COORDSNEW(:)=P(:)225:      COORDSNEW(:)=P(:)
 226:      PERCT=.FALSE.
243: !    IF (DEBUG) WRITE(MYUNIT,'(3G20.10)') P(1:3*NATOMS)227: !    IF (DEBUG) WRITE(MYUNIT,'(3G20.10)') P(1:3*NATOMS)
244: !     IF (DEBUG) THEN !this debug printing does not make sense anymore, we already reordered the coordinates228: !     IF (DEBUG) THEN !this debug printing does not make sense anymore, we already reordered the coordinates
245: !        DO J1=1,NSITES229: !        DO J1=1,NSITES
246: !           WRITE(MYUNIT,'(2I6,3G20.10)') J1,NDIST1(J1),P(3*(J1-1)+1),P(3*(J1-1)+2),P(3*(J1-1)+3)230: !           WRITE(MYUNIT,'(2I6,3G20.10)') J1,NDIST1(J1),P(3*(J1-1)+1),P(3*(J1-1)+2),P(3*(J1-1)+3)
247: !        ENDDO231: !        ENDDO
248: !     ENDIF232: !     ENDIF
249:   ELSE !NUNCON1==0, so every atom is connected in the last group, or the very last atom is disconnected233:   ELSE !NUNCON1==0, so every atom is connected in the last group, or the very last atom is disconnected
250:    atomingroup(NSITESSKIP+1:NSITES)=groupindex234:    atomingroup(NSITESSKIP+1:NSITES)=groupindex
251:      IF(DEBUG) THEN235:      IF(DEBUG) THEN
252:         WRITE(MYUNIT,*) 'percgroup> finished grouping, particles in the following groups: '236:         WRITE(MYUNIT,*) 'percgroup> finished grouping, particles in the following groups: '


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0