hdiff output

r33132/benzgenrigid.f90 2017-08-07 15:31:02.309679859 +0100 r33131/benzgenrigid.f90 2017-08-07 15:31:05.273719286 +0100
  1: ! -----------------------------------------------------------------------------  1: ! dj337: Anisotropic potential for polycyclic aromatic hydrocarbons.
  2: ! dj337 
  3: ! Anisotropic potential for periodic benzene systems. Potential from: 
  4: ! T.S.Totton, A.J.Misquitta, M.Kraft, J. Chem. Theory Comput. 6 (2010) 683–695. 
  5:  
  6: ! Long-range electrostatic interactions are computed using Ewald summation.  2: ! Long-range electrostatic interactions are computed using Ewald summation.
  7:   3: ! Implemented within the GENRIGID framework.
  8: ! Implemented within the GENRIGID framework for rigid bodies. 
  9:  
 10: ! Gradients wrt lattice parameters are computed when BOXDERIV keyword used. This 
 11: ! subroutine is for triclinic systems. Use BENZGENRIGIDEWALD_ORTHO for orthorhombic 
 12: ! cells. 
 13: ! ----------------------------------------------------------------------------- 
 14:   4: 
 15:       SUBROUTINE BENZGENRIGIDEWALD(X, G, ENERGY, GTEST)  5:       SUBROUTINE BENZGENRIGIDEWALD(X, G, ENERGY, GTEST)
 16:   6: 
 17:       USE COMMONS, ONLY: NATOMS, NCARBON, RBSTLA, RHOCC0, RHOCC10, RHOCC20, &  7:       USE COMMONS, ONLY: NATOMS, NCARBON, RBSTLA, RHOCC0, RHOCC10, RHOCC20, &
 18:      &                   RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, RHOCH10, RHOC20H, &  8:      &                   RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, RHOCH10, RHOC20H, &
 19:      &                   RHOCH20, ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, &  9:      &                   RHOCH20, ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, &
 20:      &                   EWALDREALC, BOX_PARAMS, BOX_PARAMSGRAD 10:      &                   EWALDREALC, BOX_PARAMS, BOX_PARAMSGRAD
 21:  11: 
 22:       ! adapted to the genrigid framework 12:       ! dj337: PAHA adapted to the genrigid framework
 23:       USE GENRIGID, ONLY: NRIGIDBODY, ATOMRIGIDCOORDT, TRANSFORMCTORIGID, NSITEPERBODY, & 13:       USE GENRIGID, ONLY: NRIGIDBODY, ATOMRIGIDCOORDT, TRANSFORMCTORIGID, NSITEPERBODY, &
 24:      &                    MAXSITE, SITESRIGIDBODY, TRANSFORMRIGIDTOC, TRANSFORMGRAD, INVERSEMATRIX 14:      &                    MAXSITE, SITESRIGIDBODY, TRANSFORMRIGIDTOC, TRANSFORMGRAD, INVERSEMATRIX
 25:  15: 
 26:       ! use Ewald summation to compute electrostatics 16:       ! dj337: use Ewald summation to compute electrostatics
 27:       USE EWALD 17:       USE EWALD
 28:       ! fractional coordinates, triclinic cells, box derivatives 
 29:       USE CARTDIST 18:       USE CARTDIST
 30:       USE BOX_DERIVATIVES 19:       USE BOX_DERIVATIVES
 31:  20: 
 32:       IMPLICIT NONE 21:       IMPLICIT NONE
 33:  22: 
 34:       INTEGER          :: I, J, K, J1, J2, J3, J4, J5, J6, J7, J8, OFFSET, FCT(6), L, M, N, IDX 23:       INTEGER          :: I, J, K, J1, J2, J3, J4, J5, J6, J7, J8, OFFSET, FCT(6), L, M, N, IDX
 35:       INTEGER          :: NEWALDREAL(3) 24:       INTEGER          :: NEWALDREAL(3)
 36:       DOUBLE PRECISION :: X(3*NATOMS) 25:       DOUBLE PRECISION :: X(3*NATOMS)
 37:       DOUBLE PRECISION, INTENT(OUT) :: G(3*NATOMS) 26:       DOUBLE PRECISION, INTENT(OUT) :: G(3*NATOMS)
 38:       DOUBLE PRECISION :: XR(3*NATOMS), XC(3*NATOMS), G3C(3*NATOMS), G3(3*NATOMS) 27:       DOUBLE PRECISION :: XR(3*NATOMS), XC(3*NATOMS), G3C(3*NATOMS), G3(3*NATOMS)
 47:       DOUBLE PRECISION :: RHOCC, RHOHH, RHOCH, COSTA, COSTB, DMPFCT, DDMPDR, EXPFCT  36:       DOUBLE PRECISION :: RHOCC, RHOHH, RHOCH, COSTA, COSTB, DMPFCT, DDMPDR, EXPFCT 
 48:       DOUBLE PRECISION :: DRIJDPI(3), DRIJDPJ(3), DCADPI(3), DCBDPI(3), DCADPJ(3), DCBDPJ(3) 37:       DOUBLE PRECISION :: DRIJDPI(3), DRIJDPJ(3), DCADPI(3), DCBDPI(3), DCADPJ(3), DCBDPJ(3)
 49:       DOUBLE PRECISION :: H(3,3), H_grad(3,3,6), H_inverse(3,3), rrfrac(3), rssfracmin(3), rrcom(3), rcomfrac(3) 38:       DOUBLE PRECISION :: H(3,3), H_grad(3,3,6), H_inverse(3,3), rrfrac(3), rssfracmin(3), rrcom(3), rcomfrac(3)
 50:       double precision :: rrcomfrac(3), rcomfracmin(3), rssfrac(3), vol, v_fact, dv_fact(3), c(3), s(3) 39:       double precision :: rrcomfrac(3), rcomfracmin(3), rssfrac(3), vol, v_fact, dv_fact(3), c(3), s(3)
 51:       double precision :: reciplatvec(3,3), reciplatvec_grad(3,3,6) 40:       double precision :: reciplatvec(3,3), reciplatvec_grad(3,3,6)
 52:       DOUBLE PRECISION, PARAMETER :: B = 1.6485D0 41:       DOUBLE PRECISION, PARAMETER :: B = 1.6485D0
 53:       double precision, parameter :: pi = 3.141592654d0 42:       double precision, parameter :: pi = 3.141592654d0
 54:       integer, parameter          :: image_cutoff = 5 43:       integer, parameter          :: image_cutoff = 5
 55:       LOGICAL          :: GTEST, keep_angles 44:       LOGICAL          :: GTEST, keep_angles
 56:  45: 
 57:       ! check if combination of triclinic cell angles is realistic 46: 
 58:       if (boxderivt) then 47:       if (boxderivt) then
 59:          keep_angles = check_angles(box_params(4:6)) 48:          keep_angles = check_angles(box_params(4:6))
 60:          if (.not.keep_angles) then 49:          if (.not.keep_angles) then
 61:             ! reject the structure if unrealistic cell angles 50:             print *, 'rejecting1'
 62:             call reject(energy, g) 51:             call reject(energy, g)
 63:             return 52:             return
 64:          endif 53:          endif
 65:       endif 54:       endif
 66:  55: 
 67:       call build_H(H, H_grad, gtest) ! H matrix 56:       call build_H(H, H_grad, gtest)
 68:       call inversematrix(H, H_inverse) ! inverse of H matrix 57:       call inversematrix(H, H_inverse)
 69:       call get_volume(vol) ! cell volume 58:       call get_volume(vol)
 70:       call get_reciplatvec(reciplatvec,reciplatvec_grad, .false.) ! reciprocal lattice vectors 59:       call get_reciplatvec(reciplatvec,reciplatvec_grad, .false.)
 71:  60: 
 72:       ! figure out how many lattice vectors to sum over 
 73:       newaldreal(1) = floor(ewaldrealc*dsqrt(sum(reciplatvec(1,:)**2))/(2.0d0*pi) + 0.5d0) 61:       newaldreal(1) = floor(ewaldrealc*dsqrt(sum(reciplatvec(1,:)**2))/(2.0d0*pi) + 0.5d0)
 74:       newaldreal(2) = floor(ewaldrealc*dsqrt(sum(reciplatvec(2,:)**2))/(2.0d0*pi) + 0.5d0) 62:       newaldreal(2) = floor(ewaldrealc*dsqrt(sum(reciplatvec(2,:)**2))/(2.0d0*pi) + 0.5d0)
 75:       newaldreal(3) = floor(ewaldrealc*dsqrt(sum(reciplatvec(3,:)**2))/(2.0d0*pi) + 0.5d0) 63:       newaldreal(3) = floor(ewaldrealc*dsqrt(sum(reciplatvec(3,:)**2))/(2.0d0*pi) + 0.5d0)
  64:       !print *, 'newaldreal: ', newaldreal(:3)
 76:  65: 
 77:       ! reject structure if would have to sum over more than five lattice vectors 
 78:       ! NOTE: this is because cells have tendency to flatten out... it takes a long time 
 79:       ! to compute these structures and there are probably equivalent ones of a more 
 80:       ! regular shape 
 81:       if (boxderivt) then 66:       if (boxderivt) then
 82:          if (.not. all(newaldreal.le.image_cutoff)) then 67:          if (.not. all(newaldreal.le.image_cutoff)) then
  68:             print *, 'rejecting2'
 83:             call reject(energy, g) 69:             call reject(energy, g)
 84:             return 70:             return
 85:          endif 71:          endif
 86:       endif 72:       endif
 87:  73: 
 88:       ! factorials 74:       ! factorials
 89:       FCT(1) = 1; FCT(2) = 2; FCT(3) = 6; FCT(4) = 24; FCT(5) = 120; FCT(6) = 720 75:       FCT(1) = 1; FCT(2) = 2; FCT(3) = 6; FCT(4) = 24; FCT(5) = 120; FCT(6) = 720
 90:       ! initialize energy values 76:       ! initialize energy values
 91:       ! energy1 is due to short-range anisotropic interactions 77:       ! energy1 is due to short-range anisotropic interactions
 92:       ! energy2 is due to damped dispersion 78:       ! energy2 is due to damped dispersion
 93:       ! energy3 is due to long-range electrostatics (computed using Ewald) 79:       ! energy3 is due to long-range electrostatics (computed using Ewald)
 94:       ENERGY = 0.D0; ENERGY1 = 0.D0; ENERGY2 = 0.D0; ENERGY3 = 0.D0 80:       ENERGY = 0.D0; ENERGY1 = 0.D0; ENERGY2 = 0.D0; ENERGY3 = 0.D0
 95:  81: 
 96:       ! initialize gradient if GTEST true 82:       ! initialize gradient if GTEST true
 97:       IF (GTEST) G(:) = 0.D0 83:       IF (GTEST) G(:) = 0.D0
 98:       IF (GTEST) G3C(:) = 0.D0 84:       IF (GTEST) G3C(:) = 0.D0
 99:  85: 
  86:       !print *, 'boxparams benzrigid: ', box_params(1:6)
  87:       !print *, 'coords benzrigid   : ', x(1:25)
  88: 
100:       ! dj337: check if input coordinates are cartesian 89:       ! dj337: check if input coordinates are cartesian
101:       ! assumes ATOMRIGIDCOORDT is correct 90:       ! assumes ATOMRIGIDCOORDT is correct
102:       IF (ATOMRIGIDCOORDT) THEN ! if input is cartesian 91:       IF (ATOMRIGIDCOORDT) THEN ! if input is cartesian
  92:          !print *, 'converting...'
103:          ! convert to rigidbody coordinates 93:          ! convert to rigidbody coordinates
104:          XR(:) = 0.D0 94:          XR(:) = 0.D0
105:          CALL TRANSFORMCTORIGID(X, XR) 95:          CALL TRANSFORMCTORIGID(X, XR)
106:          if (boxderivt) then 96:          if (boxderivt) then
107:             call frac2cart_rb_tri(nrigidbody, xdum, xr, H) 97:             call frac2cart_rb_tri(nrigidbody, xdum, xr, H)
  98:             !if (ortho) call frac2cart_rb_ortho(xdum, xr)
108:             x(:) = xdum(:) 99:             x(:) = xdum(:)
109:          else100:          else
110:             x(:) = xr(:)101:             x(:) = xr(:)
111:          endif102:          endif
112:       ENDIF103:       ENDIF
113: 104: 
 105:       !print *, 'coords benzrigid   : '
 106:       !do j1 = 1, 2*nrigidbody
 107:       !   print *, x(3*j1-2:3*j1)
 108:       !enddo
 109: 
 110:       !call build_H(H, H_grad, gtest)
 111:       !call inversematrix(H, H_inverse)
 112:       !call get_volume(vol)
 113:       !call get_reciplatvec(reciplatvec,reciplatvec_grad, .false.)
114:       if (boxderivt) then114:       if (boxderivt) then
115:          ! compute v factor115:          ! compute v factor
116:          c(:) = dcos(box_params(4:6))116:          c(:) = dcos(box_params(4:6))
117:          s(:) = dsin(box_params(4:6))117:          s(:) = dsin(box_params(4:6))
118:          ! v_fact is factor related to volume 
119:          v_fact = dsqrt(1.0d0 - c(1)**2-c(2)**2-c(3)**2 + 2.0d0*c(1)*c(2)*c(3))118:          v_fact = dsqrt(1.0d0 - c(1)**2-c(2)**2-c(3)**2 + 2.0d0*c(1)*c(2)*c(3))
120:          dv_fact(1) = s(1)*(c(1) - c(2)*c(3))/v_fact119:          dv_fact(1) = s(1)*(c(1) - c(2)*c(3))/v_fact
121:          dv_fact(2) = s(2)*(c(2) - c(1)*c(3))/v_fact120:          dv_fact(2) = s(2)*(c(2) - c(1)*c(3))/v_fact
122:          dv_fact(3) = s(3)*(c(3) - c(1)*c(2))/v_fact121:          dv_fact(3) = s(3)*(c(3) - c(1)*c(2))/v_fact
123:       endif122:       endif
124: 123: 
125:       EWALDREALC2 = EWALDREALC**2 ! real-space cutoff124:       EWALDREALC2 = EWALDREALC**2
 125: 
 126:       !print *, 'reciplatvec: ', reciplatvec(:3,:3)
126: 127: 
127:       ! OFFSET is number of CoM coords (3*NRIGIDBODY)128:       ! OFFSET is number of CoM coords (3*NRIGIDBODY)
128:       OFFSET     = 3*NRIGIDBODY129:       OFFSET     = 3*NRIGIDBODY
129: 130: 
130:       ! Computing Cartesian coordinates for the system.  131:       ! Computing Cartesian coordinates for the system.  
131:       DO J1 = 1, NRIGIDBODY132:       DO J1 = 1, NRIGIDBODY
132: 133: 
133:          J3 = 3*J1134:          J3 = 3*J1
134:          J5 = OFFSET + J3135:          J5 = OFFSET + J3
135:          ! CoM coords for rigid body J1136:          ! CoM coords for rigid body J1
157:                ! calculate derivative wrt coordinates158:                ! calculate derivative wrt coordinates
158:                DR1(J4,:) = MATMUL(DRMI1(:,:),SITESRIGIDBODY(J2,:,J1))159:                DR1(J4,:) = MATMUL(DRMI1(:,:),SITESRIGIDBODY(J2,:,J1))
159:                DR2(J4,:) = MATMUL(DRMI2(:,:),SITESRIGIDBODY(J2,:,J1))160:                DR2(J4,:) = MATMUL(DRMI2(:,:),SITESRIGIDBODY(J2,:,J1))
160:                DR3(J4,:) = MATMUL(DRMI3(:,:),SITESRIGIDBODY(J2,:,J1))161:                DR3(J4,:) = MATMUL(DRMI3(:,:),SITESRIGIDBODY(J2,:,J1))
161: 162: 
162:                ! calculate derivative wrt local axis163:                ! calculate derivative wrt local axis
163:                DE1(J4,:) = MATMUL(DRMI1(:,:),RBSTLA(J2,:))164:                DE1(J4,:) = MATMUL(DRMI1(:,:),RBSTLA(J2,:))
164:                DE2(J4,:) = MATMUL(DRMI2(:,:),RBSTLA(J2,:))165:                DE2(J4,:) = MATMUL(DRMI2(:,:),RBSTLA(J2,:))
165:                DE3(J4,:) = MATMUL(DRMI3(:,:),RBSTLA(J2,:))166:                DE3(J4,:) = MATMUL(DRMI3(:,:),RBSTLA(J2,:))
166: 167: 
167:             ENDIF ! gtest168:             ENDIF
 169: 
 170:          ENDDO
168: 171: 
169:          ENDDO ! sites172:       ENDDO
170: 173: 
171:       ENDDO ! rigid bodies174:       !print *, 'cart coords benzrigid   : '
 175:       !do j1 = 1, natoms
 176:       !   print *, r(j1,:3)
 177:       !enddo
172: 178: 
173:       ! Now compute the actual potential.179:       ! Now compute the actual potential.
174:       ! loop over rigid bodies (A)180:       ! loop over rigid bodies (A)
175:       DO J1 = 1, NRIGIDBODY - 1181:       DO J1 = 1, NRIGIDBODY - 1
176: 182: 
177:          J3 = 3*J1183:          J3 = 3*J1
178:          J5 = OFFSET + J3184:          J5 = OFFSET + J3
179:          ! CoM coords for rigid body J1185:          ! CoM coords for rigid body J1
180:          RI(:)  = X(J3-2:J3)186:          RI(:)  = X(J3-2:J3)
181: 187: 
212:                      ! get center of mass separation vector218:                      ! get center of mass separation vector
213:                      rrcom(:) = x(j3-2:j3) - x(j4-2:j4)219:                      rrcom(:) = x(j3-2:j3) - x(j4-2:j4)
214:                      ! convert to fractional coordinates220:                      ! convert to fractional coordinates
215:                      rrcomfrac(:) = matmul(H_inverse, rrcom(:))221:                      rrcomfrac(:) = matmul(H_inverse, rrcom(:))
216:                      ! minimum image convention222:                      ! minimum image convention
217:                      rcomfracmin(1) = rrcomfrac(1) - anint(rrfrac(1))223:                      rcomfracmin(1) = rrcomfrac(1) - anint(rrfrac(1))
218:                      rcomfracmin(2) = rrcomfrac(2) - anint(rrfrac(2))224:                      rcomfracmin(2) = rrcomfrac(2) - anint(rrfrac(2))
219:                      rcomfracmin(3) = rrcomfrac(3) - anint(rrfrac(3))225:                      rcomfracmin(3) = rrcomfrac(3) - anint(rrfrac(3))
220:                   endif226:                   endif
221: 227: 
222:                   ! sum over lattice vectors 
223:                   do l = -newaldreal(1), newaldreal(1)228:                   do l = -newaldreal(1), newaldreal(1)
224:                   rssfrac(1) = rssfracmin(1) + l229:                   rssfrac(1) = rssfracmin(1) + l
225: 230: 
226:                      do m = -newaldreal(2), newaldreal(2)231:                      do m = -newaldreal(2), newaldreal(2)
227:                      rssfrac(2) = rssfracmin(2) + m232:                      rssfrac(2) = rssfracmin(2) + m
228: 233: 
229:                         do n = -newaldreal(3), newaldreal(3)234:                         do n = -newaldreal(3), newaldreal(3)
230:                         rssfrac(3) = rssfracmin(3) + n235:                         rssfrac(3) = rssfracmin(3) + n
231: 236: 
232:                         ! convert to absolute coordinates237:                         ! convert to absolute coordinates
233:                         rss(:) = matmul(H, rssfrac(:))238:                         rss(:) = matmul(H, rssfrac(:))
234: 239: 
235:                         ! get COM vector240:                         ! get COM vector
236:                         if (gtest.and.boxderivt) then241:                         if (gtest.and.boxderivt) then
237:                            rcomfrac(1) = rcomfracmin(1) + l242:                            rcomfrac(1) = rcomfracmin(1) + l
238:                            rcomfrac(2) = rcomfracmin(2) + m243:                            rcomfrac(2) = rcomfracmin(2) + m
239:                            rcomfrac(3) = rcomfracmin(3) + n244:                            rcomfrac(3) = rcomfracmin(3) + n
240:                         endif245:                         endif
241:                      246:                      
242:                         R2     = DOT_PRODUCT(RSS(:),RSS(:))247:                         R2     = DOT_PRODUCT(RSS(:),RSS(:))
 248:                         !print *, 'r2: ', r2
243:                         ! check if distance within cutoff249:                         ! check if distance within cutoff
244:                         IF (R2 < EWALDREALC2) THEN250:                         IF (R2 < EWALDREALC2) THEN
 251:                            !print *, j7, j8
 252:                            !print *, 'r   : ', rss(:3)
 253:                            !print *, 'rmin: ', rssmin(1:3)
 254:                            !print *, 'rcom: ', rcom(1:3)
245:                            ! ABSRIJ is site-site separation between I and J255:                            ! ABSRIJ is site-site separation between I and J
246:                            ABSRIJ = DSQRT(R2)256:                            ABSRIJ = DSQRT(R2)
247:                            ! NR is unit site-site vector from sites I to J257:                            ! NR is unit site-site vector from sites I to J
248:                            NR(:)  = RSS(:)/ABSRIJ258:                            NR(:)  = RSS(:)/ABSRIJ
249:                            R2     = 1.D0/R2259:                            R2     = 1.D0/R2
250:                            R6     = R2*R2*R2260:                            R6     = R2*R2*R2
 261:                            !print *, 'absrij: ', absrij
251:          262:          
252:       !     CALCULATE THE DISPERSION DAMPING FACTOR263:       !     CALCULATE THE DISPERSION DAMPING FACTOR
253:          264:          
254:                            ! initialize sum for the damping function and vertical shift265:                            ! initialize sum for the damping function and vertical shift
255:                            DMPFCT = 1.D0266:                            DMPFCT = 1.D0
256:                            DMPFCT_SHIFT = 1.D0267:                            DMPFCT_SHIFT = 1.D0
257:                            ! initialize sum for the derivative of damping function268:                            ! initialize sum for the derivative of damping function
258:                            DDMPDR = B269:                            DDMPDR = B
259:          270:          
260:                            ! calculate sums271:                            ! calculate sums
330:                               RHOCC   = RHOCC0 + RHOCC10*(COSTA + COSTB) + RHOCC20*(1.5D0*COSTA*COSTA & 341:                               RHOCC   = RHOCC0 + RHOCC10*(COSTA + COSTB) + RHOCC20*(1.5D0*COSTA*COSTA & 
331:                                       + 1.5D0*COSTB*COSTB - 1.D0)342:                                       + 1.5D0*COSTB*COSTB - 1.D0)
332:                               ! ENERGY1 is energy due to short-range anisotropic interactions343:                               ! ENERGY1 is energy due to short-range anisotropic interactions
333:                               ! calculate vertical shift for first term344:                               ! calculate vertical shift for first term
334:                               EXPFCT  = KKJ*DEXP(-ALPHACC*(ABSRIJ - RHOCC))345:                               EXPFCT  = KKJ*DEXP(-ALPHACC*(ABSRIJ - RHOCC))
335:                               VSHIFT1 = KKJ*DEXP(-ALPHACC*(EWALDREALC - RHOCC))346:                               VSHIFT1 = KKJ*DEXP(-ALPHACC*(EWALDREALC - RHOCC))
336:                               ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1347:                               ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
337:                               ! ENERGY2 is energy due to damped dispersion348:                               ! ENERGY2 is energy due to damped dispersion
338:                               ! calculate vertical shift for second term349:                               ! calculate vertical shift for second term
339:                               VSHIFT2 = DC6CC*DMPFCT_SHIFT/(EWALDREALC**6)350:                               VSHIFT2 = DC6CC*DMPFCT_SHIFT/(EWALDREALC**6)
 351:                               !print *, 'energy: ', dc6cc*dmpfct*r6
340:                               ENERGY2 = ENERGY2 - DC6CC*DMPFCT*R6 + VSHIFT2352:                               ENERGY2 = ENERGY2 - DC6CC*DMPFCT*R6 + VSHIFT2
341:          353:          
342:                               IF (GTEST) THEN354:                               IF (GTEST) THEN
343:          355:          
344:                                  ! DVDR is derivative of dispersion damping factor energy with factor 1/Rab356:                                  ! DVDR is derivative of dispersion damping factor energy with factor 1/Rab
345:                                  DVDR    = 6.D0*DC6CC*R6*R2*DMPFCT - DC6CC*R6*DDMPDR 357:                                  DVDR    = 6.D0*DC6CC*R6*R2*DMPFCT - DC6CC*R6*DDMPDR 
 358:                                  !print *, 'grad: ', dvdr
346:                                  ! FRIJ is derivative of ENERGY1 wrt r_ij with factor 1/Rab359:                                  ! FRIJ is derivative of ENERGY1 wrt r_ij with factor 1/Rab
347:                                  FRIJ(:) = ALPHACC*EXPFCT*(-NR(:) + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADR(:) &360:                                  FRIJ(:) = ALPHACC*EXPFCT*(-NR(:) + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADR(:) &
348:                                          + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDR(:))361:                                          + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDR(:))
349:                                  ! TIJ is derivative of ENERGY1 wrt pi with factor 1/Rab362:                                  ! TIJ is derivative of ENERGY1 wrt pi with factor 1/Rab
350:                                  TIJ(:)  = ALPHACC*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADPI(:) &363:                                  TIJ(:)  = ALPHACC*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADPI(:) &
351:                                          + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDPI(:))364:                                          + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDPI(:))
352:                                  ! TJI is derivative of ENERGY1 wrt pj with factor 1/Rab365:                                  ! TJI is derivative of ENERGY1 wrt pj with factor 1/Rab
353:                                  TJI(:)  = ALPHACC*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADPJ(:) &366:                                  TJI(:)  = ALPHACC*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADPJ(:) &
354:                                          + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDPJ(:)) 367:                                          + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDPJ(:)) 
355:          368:          
357:          370:          
358:                            ! calculate if I and J are both hydorgens371:                            ! calculate if I and J are both hydorgens
359:                            ELSEIF (I > NCARBON .AND. J > NCARBON) THEN372:                            ELSEIF (I > NCARBON .AND. J > NCARBON) THEN
360:          373:          
361:                               RHOHH  = RHOHH0 + RHOHH10*(COSTA + COSTB) + RHOHH20*(1.5D0*COSTA*COSTA      &374:                               RHOHH  = RHOHH0 + RHOHH10*(COSTA + COSTB) + RHOHH20*(1.5D0*COSTA*COSTA      &
362:                                      + 1.5D0*COSTB*COSTB - 1.D0) 375:                                      + 1.5D0*COSTB*COSTB - 1.D0) 
363:                               EXPFCT  = KKJ*DEXP(-ALPHAHH*(ABSRIJ - RHOHH))376:                               EXPFCT  = KKJ*DEXP(-ALPHAHH*(ABSRIJ - RHOHH))
364:                               VSHIFT1 = KKJ*DEXP(-ALPHAHH*(EWALDREALC - RHOHH))377:                               VSHIFT1 = KKJ*DEXP(-ALPHAHH*(EWALDREALC - RHOHH))
365:                               ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1378:                               ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
366:                               VSHIFT2 = DC6HH*DMPFCT_SHIFT/(EWALDREALC**6)379:                               VSHIFT2 = DC6HH*DMPFCT_SHIFT/(EWALDREALC**6)
 380:                               !print *, 'energy: ', dc6hh*dmpfct*r6
367:                               ENERGY2 = ENERGY2 - DC6HH*DMPFCT*R6 + VSHIFT2381:                               ENERGY2 = ENERGY2 - DC6HH*DMPFCT*R6 + VSHIFT2
368:          382:          
369:                               IF (GTEST) THEN383:                               IF (GTEST) THEN
370:          384:          
371:                                  DVDR    = 6.D0*DC6HH*R6*R2*DMPFCT - DC6HH*R6*DDMPDR 385:                                  DVDR    = 6.D0*DC6HH*R6*R2*DMPFCT - DC6HH*R6*DDMPDR 
 386:                                  !print *, 'grad: ', dvdr
372:                                  FRIJ(:) = ALPHAHH*EXPFCT*(-NR(:) + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADR(:) &387:                                  FRIJ(:) = ALPHAHH*EXPFCT*(-NR(:) + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADR(:) &
373:                                          + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDR(:))388:                                          + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDR(:))
374:                                  TIJ(:)  = ALPHAHH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPI(:) &389:                                  TIJ(:)  = ALPHAHH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPI(:) &
375:                                          + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPI(:))390:                                          + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPI(:))
376:                                  TJI(:)  = ALPHAHH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPJ(:) &391:                                  TJI(:)  = ALPHAHH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPJ(:) &
377:                                          + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPJ(:))392:                                          + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPJ(:))
378:          393:          
379:                               ENDIF394:                               ENDIF
380:          395:          
381:                            ! calculate if I is carbon and J is hydrogen396:                            ! calculate if I is carbon and J is hydrogen
382:                            ELSE IF (I <= NCARBON .AND. J > NCARBON) THEN 397:                            ELSE IF (I <= NCARBON .AND. J > NCARBON) THEN 
383:          398:          
384:                               RHOCH  = RHOCH0 + RHOC10H*COSTA + RHOCH10*COSTB + RHOC20H*(1.5D0*COSTA*COSTA &399:                               RHOCH  = RHOCH0 + RHOC10H*COSTA + RHOCH10*COSTB + RHOC20H*(1.5D0*COSTA*COSTA &
385:                                      - 0.5D0) + RHOCH20*(1.5D0*COSTB*COSTB - 0.5D0)400:                                      - 0.5D0) + RHOCH20*(1.5D0*COSTB*COSTB - 0.5D0)
386:                               EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))401:                               EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))
387:                               VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))402:                               VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))
388:                               ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1403:                               ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
389:                               VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)404:                               VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)
 405:                               !print *, 'energy: ', dc6ch*dmpfct*r6
390:                               ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2406:                               ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2
391:          407:          
392:                               IF (GTEST) THEN408:                               IF (GTEST) THEN
393:                            409:                            
394:                                  DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 410:                                  DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 
 411:                                  !print *, 'grad: ', dvdr
395:                                  FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADR(:) &412:                                  FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADR(:) &
396:                                          + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDR(:))413:                                          + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDR(:))
397:                                  TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPI(:) &414:                                  TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPI(:) &
398:                                          + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPI(:))415:                                          + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPI(:))
399:                                  TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPJ(:) &416:                                  TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPJ(:) &
400:                                          + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPJ(:))417:                                          + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPJ(:))
401:          418:          
402:                               ENDIF419:                               ENDIF
403:          420:          
404:                            ELSE !IF(I > NCARBON .AND. J <= NCARBON) THEN421:                            ELSE !IF(I > NCARBON .AND. J <= NCARBON) THEN
405:          422:          
406:                               RHOCH  = RHOCH0 + RHOCH10*COSTA + RHOC10H*COSTB + RHOCH20*(1.5D0*COSTA*COSTA &423:                               RHOCH  = RHOCH0 + RHOCH10*COSTA + RHOC10H*COSTB + RHOCH20*(1.5D0*COSTA*COSTA &
407:                                      - 0.5D0) + RHOC20H*(1.5D0*COSTB*COSTB - 0.5D0)424:                                      - 0.5D0) + RHOC20H*(1.5D0*COSTB*COSTB - 0.5D0)
408:                               EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))425:                               EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))
409:                               VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))426:                               VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))
410:                               ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1427:                               ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
411:                               VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)428:                               VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)
 429:                               !print *, 'energy: ', dc6ch*dmpfct*r6
412:                               ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2430:                               ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2
413:          431:          
414:                               IF (GTEST) THEN432:                               IF (GTEST) THEN
415:          433:          
416:                                  DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 434:                                  DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 
 435:                                  !print *, 'grad: ', dvdr
417:                                  FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADR(:) &436:                                  FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADR(:) &
418:                                          + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDR(:))437:                                          + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDR(:))
419:                                  TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPI(:) &438:                                  TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPI(:) &
420:                                          + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPI(:))439:                                          + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPI(:))
421:                                  TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPJ(:) &440:                                  TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPJ(:) &
422:                                          + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPJ(:))441:                                          + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPJ(:))
423:          442:          
424:                               ENDIF443:                               ENDIF
425:          444:          
426:                            ENDIF445:                            ENDIF
427:          446:          
428:                            IF (GTEST) THEN447:                            IF (GTEST) THEN
429:          448:          
430:                               ! total gradient wrt CoM coords for rigid body J1449:                               ! total gradient wrt CoM coords for rigid body J1
431:                               G(J3-2:J3) = G(J3-2:J3) + DVDR*RSS(:) + FRIJ(:)450:                               G(J3-2:J3) = G(J3-2:J3) + DVDR*RSS(:) + FRIJ(:)
 451:                               !g(j3-2:j3) = g(j3-2:j3) + frij(:)
 452: 
 453:                               !box_paramsgrad(1:3) = box_paramsgrad(1:3) + (dvdr*rss(1:3)+frij(1:3))*rcom(1:3)/box_params(1:3)
432:                               ! total gradient wrt CoM coords for rigid body J2454:                               ! total gradient wrt CoM coords for rigid body J2
433:                               G(J4-2:J4) = G(J4-2:J4) - DVDR*RSS(:) - FRIJ(:)455:                               G(J4-2:J4) = G(J4-2:J4) - DVDR*RSS(:) - FRIJ(:)
 456:                               !g(j4-2:j4) = g(j4-2:j4) - frij(:)
434: 457: 
435:                               ! total gradient wrt AA coords for rigid body J1458:                               ! total gradient wrt AA coords for rigid body J1
436:                               G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPI(:) + TIJ(:)459:                               G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPI(:) + TIJ(:)
 460:                               !g(j5-2:j5) = g(j5-2:j5) + tij(:)
437:                               ! total gradient wrt AA coords for rigid body J2461:                               ! total gradient wrt AA coords for rigid body J2
438:                               G(J6-2:J6) = G(J6-2:J6) + DVDR*DRIJDPJ(:) + TJI(:)462:                               G(J6-2:J6) = G(J6-2:J6) + DVDR*DRIJDPJ(:) + TJI(:)
 463:                               !g(j6-2:j6) = g(j6-2:j6) + tji(:)
 464: 
 465:                               ! TODO: this is if orientation of rigid body depends on cell parameters
 466:                               !H_mat(:,:) = (-1.0d0/(vol*v_fact))*dv_fact(1)*H(:,:) + (1.0d0/vol)*H_grad(:,:,4)
 467:                               !rmatrix(:) = matmul(H_inverse, vol*(rss(:)-rcom(:)))
 468:                               !box_paramsgrad(4) = box_paramsgrad(4) + dvdr*dot_product(rss(1:3), matmul(H_mat, rmatrix(:)))
439: 469: 
440:                               ! gradients wrt cell parameters470:                               !box_paramsgrad(1) = box_paramsgrad(1) + dvdr*dot_product(rss(1:3), matmul(H_mat/vol, rmatrix(:)))
441:                               if (boxderivt) then471:                               if (boxderivt) then
442:                                  do idx = 1, 6472:                                  do idx = 1, 6
443:                                     box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product((dvdr*rss(1:3) + frij(1:3)), matmul(H_grad(:,:,idx), rcomfrac(:)))473:                                     box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product((dvdr*rss(1:3) + frij(1:3)), matmul(H_grad(:,:,idx), rcomfrac(:)))
444:                                  enddo474:                                  enddo
445:                               endif ! box derivatives475:                               endif
446: 476: 
447:                            ENDIF ! gtest477:                            ENDIF
448:       478:       
449:                         ENDIF ! within cutoff479:                         ENDIF
450: 480: 
451:                      enddo ! n481:                      enddo
452:                   enddo ! m482:                   enddo
453:                enddo ! l483:                enddo
454: 484: 
455:                ENDDO ! sites j485:                ENDDO
456: 486: 
457:             ENDDO ! rigid bodies J487:             ENDDO
458:  488:  
459:          ENDDO ! sites i489:          ENDDO
460: 490: 
461:       ENDDO ! rigid bodies I491:       ENDDO
462: 492: 
463: ! INCLUDE CONTRIBUTION OF RIGID BODY WITH PERIODIC IMAGE OF ITSELF493: ! INCLUDE CONTRIBUTION OF RIGID BODY WITH PERIODIC IMAGE OF ITSELF
464: 494: 
465:       ! loop over rigidbodies495:       ! loop over rigidbodies
466:       do j1 = 1, nrigidbody496:       do j1 = 1, nrigidbody
467:          j3 = 3*j1497:          j3 = 3*j1
468:          j5 = offset + j3498:          j5 = offset + j3
469:          ri(:) = x(j3-2:j3)499:          ri(:) = x(j3-2:j3)
470: 500: 
471:          ! loop over sites i501:          ! loop over sites i
472:          do i = 1, nsiteperbody(j1)502:          do i = 1, nsiteperbody(j1)
473:             j7 = maxsite*(j1-1) + i503:             j7 = maxsite*(j1-1) + i
474:             ei(:) = e(j7,:)504:             ei(:) = e(j7,:)
475: 505: 
476:             ! loop over sites j506:             ! loop over sites j
477:             do j = 1, nsiteperbody(j1)507:             do j = 1, nsiteperbody(j1)
478:                j8 = maxsite*(j1-1) + j508:                j8 = maxsite*(j1-1) + j
479:                ej(:) = e(j8,:)509:                ej(:) = e(j8,:)
480: 510: 
481:                ! site-site separation vector511:                !print *, 'new!!', j1, i, j
482:                rr(:) = r(j7,:) - r(j8,:)512:                rr(:) = r(j7,:) - r(j8,:)
 513:                !print *, j7, r(j7,:3)
 514:                !print *, j8, r(j8,:3)
 515:                !print *, 'rr: ', rr(:3)
483:                ! convert to fractional516:                ! convert to fractional
484:                rrfrac(:) = matmul(H_inverse, rr(:))517:                rrfrac(:) = matmul(H_inverse, rr(:))
 518:                !print *, 'rrfrac: ', rrfrac(:3)
485: 519: 
486:                ! sum over lattice vectors520:                ! loop over boxes
487:                do l = -newaldreal(1), newaldreal(1)521:                do l = -newaldreal(1), newaldreal(1)
488:                   do m = -newaldreal(2), newaldreal(2)522:                   do m = -newaldreal(2), newaldreal(2)
489:                      do n = -newaldreal(3), newaldreal(3)523:                      do n = -newaldreal(3), newaldreal(3)
490: 524: 
491:                      ! if not on same rigid body 
492:                      if (.not.(l.eq.0.and.m.eq.0.and.n.eq.0)) then525:                      if (.not.(l.eq.0.and.m.eq.0.and.n.eq.0)) then
493: 526: 
494:                         rssfrac(1) = rrfrac(1) + l527:                         rssfrac(1) = rrfrac(1) + l
495:                         rssfrac(2) = rrfrac(2) + m528:                         rssfrac(2) = rrfrac(2) + m
496:                         rssfrac(3) = rrfrac(3) + n529:                         rssfrac(3) = rrfrac(3) + n
497:                         ! convert back to absolute530:                         !print *, l, m, n
 531:                         !print *, 'rssfrac: ', rssfrac(:3)
498:                         rss(:) = matmul(H, rssfrac(:))532:                         rss(:) = matmul(H, rssfrac(:))
499: 533: 
500:                         ! get vector between COM 
501:                         if (gtest.and.boxderivt) then534:                         if (gtest.and.boxderivt) then
502:                            rcomfrac(1) = l535:                            rcomfrac(1) = l
503:                            rcomfrac(2) = m536:                            rcomfrac(2) = m
504:                            rcomfrac(3) = n537:                            rcomfrac(3) = n
505:                         endif538:                         endif
506: 539: 
507:                         r2 = dot_product(rss(:), rss(:))540:                         r2 = dot_product(rss(:), rss(:))
 541:                         !print *, 'rssmin2: ', rss(1:3)
508:                         if (r2 < ewaldrealc2) then542:                         if (r2 < ewaldrealc2) then
509: 543: 
510:                         ! absolute site-site distance544:                         !print *, 'rssmin3: ', rss(1:3)
 545:                         !print *, 'r2    : ', r2
 546:                         !print *, 'rcom  : ', rcom(1:3)
511:                         absrij = dsqrt(r2)547:                         absrij = dsqrt(r2)
 548:                         !print *, j7, j8
 549:                         !print *, absrij
512:                         nr(:) = rss(:)/absrij550:                         nr(:) = rss(:)/absrij
513:                         r2 = 1.d0/r2551:                         r2 = 1.d0/r2
514:                         r6 = r2*r2*r2552:                         r6 = r2*r2*r2
515: 553: 
516:                         ! CALCULATE DISPERSION DAMPING FACTOR554:                         ! CALCULATE DISPERSION DAMPING FACTOR
517: 555: 
518:                         ! initialize sum for the damping function and vertical shift556:                         ! initialize sum for the damping function and vertical shift
519:                         DMPFCT = 1.D0557:                         DMPFCT = 1.D0
520:                         DMPFCT_SHIFT = 1.D0558:                         DMPFCT_SHIFT = 1.D0
521:                         ! initialize sum for the derivative of damping function559:                         ! initialize sum for the derivative of damping function
593:                            RHOCC   = RHOCC0 + RHOCC10*(COSTA + COSTB) + RHOCC20*(1.5D0*COSTA*COSTA & 631:                            RHOCC   = RHOCC0 + RHOCC10*(COSTA + COSTB) + RHOCC20*(1.5D0*COSTA*COSTA & 
594:                                    + 1.5D0*COSTB*COSTB - 1.D0)632:                                    + 1.5D0*COSTB*COSTB - 1.D0)
595:                            ! ENERGY1 is energy due to short-range anisotropic interactions633:                            ! ENERGY1 is energy due to short-range anisotropic interactions
596:                            ! calculate vertical shift for first term634:                            ! calculate vertical shift for first term
597:                            EXPFCT  = KKJ*DEXP(-ALPHACC*(ABSRIJ - RHOCC))635:                            EXPFCT  = KKJ*DEXP(-ALPHACC*(ABSRIJ - RHOCC))
598:                            VSHIFT1 = KKJ*DEXP(-ALPHACC*(EWALDREALC - RHOCC))636:                            VSHIFT1 = KKJ*DEXP(-ALPHACC*(EWALDREALC - RHOCC))
599:                            ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1637:                            ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
600:                            ! ENERGY2 is energy due to damped dispersion638:                            ! ENERGY2 is energy due to damped dispersion
601:                            ! calculate vertical shift for second term639:                            ! calculate vertical shift for second term
602:                            VSHIFT2 = DC6CC*DMPFCT_SHIFT/(EWALDREALC**6)640:                            VSHIFT2 = DC6CC*DMPFCT_SHIFT/(EWALDREALC**6)
 641:                            !print *, 'energy: ', dc6cc*dmpfct*r6
603:                            ENERGY2 = ENERGY2 - DC6CC*DMPFCT*R6 + VSHIFT2642:                            ENERGY2 = ENERGY2 - DC6CC*DMPFCT*R6 + VSHIFT2
 643:                            !print *, 'vshift2     : ', vshift2
 644:                            !print *, 'contribution: ', -dc6cc*dmpfct*r6
 645:                            !print *, 'energy2     : ', energy2
604: 646: 
605:                            IF (GTEST) THEN647:                            IF (GTEST) THEN
606: 648: 
607:                               ! DVDR is derivative of dispersion damping factor energy with factor 1/Rab649:                               ! DVDR is derivative of dispersion damping factor energy with factor 1/Rab
608:                               DVDR    = 6.D0*DC6CC*R6*R2*DMPFCT - DC6CC*R6*DDMPDR 650:                               DVDR    = 6.D0*DC6CC*R6*R2*DMPFCT - DC6CC*R6*DDMPDR 
609:                               ! FRIJ is derivative of ENERGY1 wrt r_ij with factor 1/Rab651:                               !print *, 'dvdr        : ', dvdr
 652:                            !   !print *, 'grad: ', dvdr
 653:                            !   ! FRIJ is derivative of ENERGY1 wrt r_ij with factor 1/Rab
610:                               FRIJ(:) = ALPHACC*EXPFCT*(-NR(:) + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADR(:) &654:                               FRIJ(:) = ALPHACC*EXPFCT*(-NR(:) + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADR(:) &
611:                                       + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDR(:))655:                                       + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDR(:))
612:                               ! TIJ is derivative of ENERGY1 wrt pi with factor 1/Rab656:                               ! TIJ is derivative of ENERGY1 wrt pi with factor 1/Rab
613:                               TIJ(:)  = ALPHACC*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADPI(:) &657:                               TIJ(:)  = ALPHACC*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADPI(:) &
614:                                       + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDPI(:))658:                                       + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDPI(:))
615:                               ! TJI is derivative of ENERGY1 wrt pj with factor 1/Rab659:                            !   ! TJI is derivative of ENERGY1 wrt pj with factor 1/Rab
616:                               TJI(:)  = ALPHACC*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADPJ(:) &660:                               TJI(:)  = ALPHACC*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADPJ(:) &
617:                                       + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDPJ(:)) 661:                                       + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDPJ(:)) 
618: 662: 
619:                            ENDIF663:                            ENDIF
620: 664: 
621:                         ! calculate if I and J are both hydorgens665:                         ! calculate if I and J are both hydorgens
622:                         ELSEIF (I > NCARBON .AND. J > NCARBON) THEN666:                         ELSEIF (I > NCARBON .AND. J > NCARBON) THEN
623: 667: 
624:                            RHOHH  = RHOHH0 + RHOHH10*(COSTA + COSTB) + RHOHH20*(1.5D0*COSTA*COSTA      &668:                            RHOHH  = RHOHH0 + RHOHH10*(COSTA + COSTB) + RHOHH20*(1.5D0*COSTA*COSTA      &
625:                                   + 1.5D0*COSTB*COSTB - 1.D0)669:                                   + 1.5D0*COSTB*COSTB - 1.D0)
626:                            EXPFCT  = KKJ*DEXP(-ALPHAHH*(ABSRIJ - RHOHH))670:                            EXPFCT  = KKJ*DEXP(-ALPHAHH*(ABSRIJ - RHOHH))
627:                            VSHIFT1 = KKJ*DEXP(-ALPHAHH*(EWALDREALC - RHOHH))671:                            VSHIFT1 = KKJ*DEXP(-ALPHAHH*(EWALDREALC - RHOHH))
628:                            ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1672:                            ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
629:                            VSHIFT2 = DC6HH*DMPFCT_SHIFT/(EWALDREALC**6)673:                            VSHIFT2 = DC6HH*DMPFCT_SHIFT/(EWALDREALC**6)
 674:                            !print *, 'energy: ', dc6hh*dmpfct*r6
630:                            ENERGY2 = ENERGY2 - DC6HH*DMPFCT*R6 + VSHIFT2675:                            ENERGY2 = ENERGY2 - DC6HH*DMPFCT*R6 + VSHIFT2
631: 676: 
632:                            IF (GTEST) THEN677:                            IF (GTEST) THEN
633: 678: 
634:                               DVDR    = 6.D0*DC6HH*R6*R2*DMPFCT - DC6HH*R6*DDMPDR 679:                               DVDR    = 6.D0*DC6HH*R6*R2*DMPFCT - DC6HH*R6*DDMPDR 
 680:                            !   !print *, 'grad: ', dvdr
635:                               FRIJ(:) = ALPHAHH*EXPFCT*(-NR(:) + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADR(:) &681:                               FRIJ(:) = ALPHAHH*EXPFCT*(-NR(:) + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADR(:) &
636:                                       + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDR(:))682:                                       + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDR(:))
637:                               TIJ(:)  = ALPHAHH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPI(:) &683:                               TIJ(:)  = ALPHAHH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPI(:) &
638:                                       + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPI(:))684:                                       + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPI(:))
639:                               TJI(:)  = ALPHAHH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPJ(:) &685:                               TJI(:)  = ALPHAHH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPJ(:) &
640:                                       + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPJ(:))686:                                       + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPJ(:))
641: 687: 
642:                            ENDIF688:                            ENDIF
643: 689: 
644:                         ! calculate if I is carbon and J is hydrogen690:                         ! calculate if I is carbon and J is hydrogen
645:                         ELSE IF (I <= NCARBON .AND. J > NCARBON) THEN 691:                         ELSE IF (I <= NCARBON .AND. J > NCARBON) THEN 
646: 692: 
647:                            RHOCH  = RHOCH0 + RHOC10H*COSTA + RHOCH10*COSTB + RHOC20H*(1.5D0*COSTA*COSTA &693:                            RHOCH  = RHOCH0 + RHOC10H*COSTA + RHOCH10*COSTB + RHOC20H*(1.5D0*COSTA*COSTA &
648:                                   - 0.5D0) + RHOCH20*(1.5D0*COSTB*COSTB - 0.5D0)694:                                   - 0.5D0) + RHOCH20*(1.5D0*COSTB*COSTB - 0.5D0)
649:                            EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))695:                            EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))
650:                            VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))696:                            VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))
651:                            ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1697:                            ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
652:                            VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)698:                            VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)
 699:                            !print *, 'energy: ', dc6ch*dmpfct*r6
653:                            ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2700:                            ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2
654: 701: 
655:                            IF (GTEST) THEN702:                            IF (GTEST) THEN
656: 703: 
657:                               DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 704:                               DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 
 705:                            !   !print *, 'grad: ', dvdr
658:                               FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADR(:) &706:                               FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADR(:) &
659:                                       + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDR(:))707:                                       + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDR(:))
660:                               TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPI(:) &708:                               TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPI(:) &
661:                                       + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPI(:))709:                                       + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPI(:))
662:                               TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPJ(:) &710:                               TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPJ(:) &
663:                                       + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPJ(:))711:                                       + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPJ(:))
664: 712: 
665:                            ENDIF713:                            ENDIF
666: 714: 
667:                         ELSE !IF(I > NCARBON .AND. J <= NCARBON) THEN715:                         ELSE !IF(I > NCARBON .AND. J <= NCARBON) THEN
668: 716: 
669:                            RHOCH  = RHOCH0 + RHOCH10*COSTA + RHOC10H*COSTB + RHOCH20*(1.5D0*COSTA*COSTA &717:                            RHOCH  = RHOCH0 + RHOCH10*COSTA + RHOC10H*COSTB + RHOCH20*(1.5D0*COSTA*COSTA &
670:                                   - 0.5D0) + RHOC20H*(1.5D0*COSTB*COSTB - 0.5D0)718:                                   - 0.5D0) + RHOC20H*(1.5D0*COSTB*COSTB - 0.5D0)
671:                            EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))719:                            EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))
672:                            VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))720:                            VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))
673:                            ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1721:                            ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
674:                            VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)722:                            VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)
 723:                            !print *, 'energy: ', dc6ch*dmpfct*r6
675:                            ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2724:                            ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2
676: 725: 
677:                            IF (GTEST) THEN726:                            IF (GTEST) THEN
678: 727: 
679:                               DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 728:                               DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 
 729:                            !   !print *, 'grad: ', dvdr
680:                               FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADR(:) &730:                               FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADR(:) &
681:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDR(:))731:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDR(:))
682:                               TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPI(:) &732:                               TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPI(:) &
683:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPI(:))733:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPI(:))
684:                               TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPJ(:) &734:                               TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPJ(:) &
685:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPJ(:))735:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPJ(:))
686: 736: 
687:                            ENDIF737:                            ENDIF
688: 738: 
689:                         ENDIF739:                         ENDIF
690: 740: 
691: 741: 
692:                         IF (GTEST) THEN742:                         IF (GTEST) THEN
693: 743: 
694:                            ! total gradient wrt AA coords for rigid body J1744:                            ! total gradient wrt AA coords for rigid body J1
695:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPI(:) + TIJ(:)745:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPI(:) + TIJ(:)
 746:                            !g(j5-2:j5) = g(j5-2:j5) + tij(:)
696:                            ! total gradient wrt AA coords for rigid body J2747:                            ! total gradient wrt AA coords for rigid body J2
697:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPJ(:) + TJI(:)748:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPJ(:) + TJI(:)
 749:                            !g(j5-2:j5) = g(j5-2:j5) + tji(:)
 750:                            !print *, 'dispersion:'
 751:                            !print *, dvdr*drijdpi(:3)
 752:                            !print *, dvdr*drijdpj(:3)
 753:                            !print *, 'energy2: ', energy2
 754:                            !print *, 'anisotropic:'
 755:                            !print *, tij(:3)
 756:                            !print *, tji(:3)
698: 757: 
699:                            ! gradients wrt cell parameters 
700:                            if (boxderivt) then758:                            if (boxderivt) then
701:                               do idx = 1, 6759:                               do idx = 1, 6
702:                                     box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product((dvdr*rss(1:3) + frij(1:3)), matmul(H_grad(:,:,idx), rcomfrac(:)))760:                                     box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product((dvdr*rss(1:3) + frij(1:3)), matmul(H_grad(:,:,idx), rcomfrac(:)))
 761:                                     !box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product(dvdr*rss(1:3), matmul(H_grad(:,:,idx), rcomfrac(:)))
703:                               enddo762:                               enddo
704:                            endif ! box derivatives763:                            endif
705: 764: 
706:                         ENDIF ! gtest765:                         ENDIF
707:                         endif ! central box766:                         endif
708:                     endif ! within cutoff767:                     endif
709:                   enddo ! n768:                   enddo
710:                enddo ! m769:                enddo
711:             enddo ! l770:             enddo
712:             enddo ! sites j771:             enddo
713:          enddo ! sites i772:          enddo
714:       enddo ! rigid bodies773:       enddo
715: 774: 
716:       ! convert to cartesian coordinates775:       ! convert to cartesian coordinates
717:       XC(:) = 0.D0776:       XC(:) = 0.D0
718:       if (boxderivt) then777:       if (boxderivt) then
719:          xdum(:) = x(:)778:          xdum(:) = x(:)
720:          call cart2frac_rb_tri(nrigidbody, xdum, x, H_inverse)779:          call cart2frac_rb_tri(nrigidbody, xdum, x, H_inverse)
721:       endif780:       endif
722:       CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XC, X)781:       CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XC, X)
723:       ! restore cartesian rigid body coordinates782:       ! restore cartesian rigid body coordinates
724:       if (boxderivt) x(:) = xdum(:)783:       if (boxderivt) x(:) = xdum(:)
725: 784: 
726:       ! ENERGY3 and G3 are energy and gradient due to electrostatics785: !      !!! ENERGY3 and G3 are energy and gradient due to electrostatics
727:       ! computed using Ewald summation786: !      !!! computed using Ewald summation
728:       CALL EWALDSUM(1, XC, G3C, ENERGY3, GTEST)787:       CALL EWALDSUM(1, XC, G3C, ENERGY3, GTEST)
729: 788: !
730:       ! convert Ewald contribution of gradient to rigidbody coordinates789: !      !!! convert Ewald contribution of gradient to rigidbody coordinates
731:       IF (GTEST) G3(:) = 0.D0790:       IF (GTEST) G3(:) = 0.D0
732:       CALL TRANSFORMGRAD(G3C, X, G3)791:       CALL TRANSFORMGRAD(G3C, X, G3)
 792: !!      !print *, 'energy2: ', energy2
 793: 
 794:       !energy = (energy1+energy2)*2625.499d0
 795:       !if (gtest) g(:) = g(:)*2625.499d0
 796:       !energy = energy2*2625.499d0
 797:       !energy = (energy1+energy3)*2625.499d0
 798:       !if (gtest) g(:) = g3(:)*2625.499d0
 799:       !ENERGY = (ENERGY1 + ENERGY2 + ENERGY3)*2625.499D0 
 800:       !IF (GTEST) G(:) = (G(:) + G3(:))*2625.499D0
 801:       !if (gtest) box_paramsgrad(1:6) = box_paramsgrad(1:6)*2625.499D0
 802:       !print *, 'box_paramsgrad: ', box_paramsgrad(1:3)
 803:       !print *, 'energy: ', energy2
733: 804: 
734:       ! dj337: if input was cartesian, convert back to cartesian805:       ! dj337: if input was cartesian, convert back to cartesian
735:       ! assumes ATOMRIGIDCOORDT is correct806:       ! assumes ATOMRIGIDCOORDT is correct
736:       IF (ATOMRIGIDCOORDT) THEN807:       IF (ATOMRIGIDCOORDT) THEN
737: 808: 
738:          ! convert to cartesian coordinates809:          ! convert to cartesian coordinates
739:          if (boxderivt) then810:          if (boxderivt) then
 811:             !print *, 'coords 708: ', x(1:3*natoms)
740:             xdum(:) = x(:)812:             xdum(:) = x(:)
741:             call cart2frac_rb_tri(nrigidbody, xdum, x, H_inverse)813:             call cart2frac_rb_tri(nrigidbody, xdum, x, H_inverse)
 814:             !if (ortho) call cart2frac_rb_ortho(xdum, x)
 815:             !print *, 'coords 711: ', x(1:3*natoms)
742:          endif816:          endif
743:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XR, X)817:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XR, X)
744:          X(:) = XR(:)818:          X(:) = XR(:)
745:       ENDIF819:       ENDIF
746: 820: 
747:       ! add WCA-style repulsion to keep cell volume away from zero 
748:       if (boxderivt) call constrain_volume(v_fact, dv_fact, energy1, box_paramsgrad(4:6), gtest)821:       if (boxderivt) call constrain_volume(v_fact, dv_fact, energy1, box_paramsgrad(4:6), gtest)
749: 822: 
750:       ! sum energies / gradients and convert to kJ/mol 
751:       ENERGY = (ENERGY1 + ENERGY2 + ENERGY3)*2625.499D0823:       ENERGY = (ENERGY1 + ENERGY2 + ENERGY3)*2625.499D0
752:       IF (GTEST) G(:) = (G(:) + G3(:))*2625.499D0824:       IF (GTEST) G(:) = (G(:) + G3(:))*2625.499D0
753:       if (gtest) box_paramsgrad(1:6) = box_paramsgrad(1:6)*2625.499D0825:       if (gtest) box_paramsgrad(1:6) = box_paramsgrad(1:6)*2625.499D0
 826:       !print *, 'coords benzrigid    : ', x(1:3*natoms)
 827:       !print *, 'box params benzrigid: ', box_params(1:6)
 828: 
 829:       !print *, 'energy: ', energy!1, energy2, energy3
 830:       !print *, 'grad: ', g(1)
 831:       !!print *, 'grad  : '
 832:       !do j1 = 1, 2*nrigidbody
 833:       !   print *, x(3*j1-2:3*j1)
 834:       !enddo
 835: 
 836:       !print *, 'contri: ', energy1, energy2, energy3
 837:       !print *, 'energy: ', energy !, box_params(1:6)
 838:       !print *, 'g1    : ', g(1:3)
 839:       !stop
754: 840: 
755:       END SUBROUTINE BENZGENRIGIDEWALD841:       END SUBROUTINE BENZGENRIGIDEWALD
756: 842: 
757: !     ----------------------------------------------------------------------------------------------843: !     ----------------------------------------------------------------------------------------------
758: !844: !
759: !      SUBROUTINE DEFPAHARIGID()845: !      SUBROUTINE DEFPAHARIGID()
760: !846: !
761: !      USE COMMONS, ONLY: RHOCC0, RHOCC10, RHOCC20,  RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, RHOCH10, RHOC20H, RHOCH20, &847: !      USE COMMONS, ONLY: RHOCC0, RHOCC10, RHOCC20,  RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, RHOCH10, RHOC20H, RHOCH20, &
762: !                         ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, CCKJ848: !                         ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, CCKJ
763: !849: !


r33132/benzgenrigid_ortho.f90 2017-08-07 15:31:02.557683151 +0100 r33131/benzgenrigid_ortho.f90 2017-08-07 15:31:05.493722209 +0100
  1: ! -----------------------------------------------------------------------------  1: ! dj337: Anisotropic potential for polycyclic aromatic hydrocarbons.
  2: ! dj337  2: ! Long-range electrostatic interactions are computed using Ewald summation.
  3: ! Anisotropic potential for periodic benzene systems. This subroutine is for  3: ! Implemented within the GENRIGID framework.
  4: ! orthorhombic cells. See BENZGENRIGIDEWALD for more info. 
  5: ! ----------------------------------------------------------------------------- 
  6:   4: 
  7:       SUBROUTINE BENZGENRIGIDEWALD_ORTHO(X, G, ENERGY, GTEST)  5:       SUBROUTINE BENZGENRIGIDEWALD_ORTHO(X, G, ENERGY, GTEST)
  8:   6: 
  9:       USE COMMONS, ONLY: NATOMS, NCARBON, RBSTLA, RHOCC0, RHOCC10, RHOCC20, &  7:       USE COMMONS, ONLY: NATOMS, NCARBON, RBSTLA, RHOCC0, RHOCC10, RHOCC20, &
 10:      &                   RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, RHOCH10, RHOC20H, &  8:      &                   RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, RHOCH10, RHOC20H, &
 11:      &                   RHOCH20, ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, &  9:      &                   RHOCH20, ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, &
 12:      &                   EWALDREALC, BOX_PARAMS, BOX_PARAMSGRAD 10:      &                   EWALDREALC, BOX_PARAMS, BOX_PARAMSGRAD
 13:  11: 
 14:       ! adapted to the genrigid framework 12:       ! dj337: PAHA adapted to the genrigid framework
 15:       USE GENRIGID, ONLY: NRIGIDBODY, ATOMRIGIDCOORDT, TRANSFORMCTORIGID, NSITEPERBODY, & 13:       USE GENRIGID, ONLY: NRIGIDBODY, ATOMRIGIDCOORDT, TRANSFORMCTORIGID, NSITEPERBODY, &
 16:      &                    MAXSITE, SITESRIGIDBODY, TRANSFORMRIGIDTOC, TRANSFORMGRAD 14:      &                    MAXSITE, SITESRIGIDBODY, TRANSFORMRIGIDTOC, TRANSFORMGRAD
 17:  15: 
 18:       ! use Ewald summation to compute electrostatics 16:       ! dj337: use Ewald summation to compute electrostatics
 19:       USE EWALD 17:       USE EWALD
 20:       USE CARTDIST 18:       USE CARTDIST
 21:       USE BOX_DERIVATIVES 19:       USE BOX_DERIVATIVES
 22:  20: 
 23:       IMPLICIT NONE 21:       IMPLICIT NONE
 24:  22: 
 25:       INTEGER          :: I, J, K, J1, J2, J3, J4, J5, J6, J7, J8, OFFSET, FCT(6), L, M, N 23:       INTEGER          :: I, J, K, J1, J2, J3, J4, J5, J6, J7, J8, OFFSET, FCT(6), L, M, N
 26:       INTEGER          :: NEWALDREAL(3) 24:       INTEGER          :: NEWALDREAL(3)
 27:       DOUBLE PRECISION :: X(3*NATOMS) 25:       DOUBLE PRECISION :: X(3*NATOMS)
 28:       DOUBLE PRECISION, INTENT(OUT) :: G(3*NATOMS) 26:       DOUBLE PRECISION, INTENT(OUT) :: G(3*NATOMS)
 34:       DOUBLE PRECISION :: R(MAXSITE*NRIGIDBODY,3), E(3*MAXSITE*NRIGIDBODY,3), xdum(3*natoms), rssmin(3) 32:       DOUBLE PRECISION :: R(MAXSITE*NRIGIDBODY,3), E(3*MAXSITE*NRIGIDBODY,3), xdum(3*natoms), rssmin(3)
 35:       DOUBLE PRECISION :: DR1(MAXSITE*NRIGIDBODY,3), DR2(MAXSITE*NRIGIDBODY,3), DR3(MAXSITE*NRIGIDBODY,3) 33:       DOUBLE PRECISION :: DR1(MAXSITE*NRIGIDBODY,3), DR2(MAXSITE*NRIGIDBODY,3), DR3(MAXSITE*NRIGIDBODY,3)
 36:       DOUBLE PRECISION :: DE1(3*MAXSITE*NRIGIDBODY,3), DE2(3*MAXSITE*NRIGIDBODY,3), DE3(3*MAXSITE*NRIGIDBODY,3) 34:       DOUBLE PRECISION :: DE1(3*MAXSITE*NRIGIDBODY,3), DE2(3*MAXSITE*NRIGIDBODY,3), DE3(3*MAXSITE*NRIGIDBODY,3)
 37:       DOUBLE PRECISION :: RMI(3,3), DRMI1(3,3), DRMI2(3,3), DRMI3(3,3), DCADR(3), DCBDR(3) 35:       DOUBLE PRECISION :: RMI(3,3), DRMI1(3,3), DRMI2(3,3), DRMI3(3,3), DCADR(3), DCBDR(3)
 38:       DOUBLE PRECISION :: RHOCC, RHOHH, RHOCH, COSTA, COSTB, DMPFCT, DDMPDR, EXPFCT, rcom(3), RRCOMMIN(3) 36:       DOUBLE PRECISION :: RHOCC, RHOHH, RHOCH, COSTA, COSTB, DMPFCT, DDMPDR, EXPFCT, rcom(3), RRCOMMIN(3)
 39:       DOUBLE PRECISION :: DRIJDPI(3), DRIJDPJ(3), DCADPI(3), DCBDPI(3), DCADPJ(3), DCBDPJ(3), rrcom(3) 37:       DOUBLE PRECISION :: DRIJDPI(3), DRIJDPJ(3), DCADPI(3), DCBDPI(3), DCADPJ(3), DCBDPJ(3), rrcom(3)
 40:       DOUBLE PRECISION, PARAMETER :: B = 1.6485D0 38:       DOUBLE PRECISION, PARAMETER :: B = 1.6485D0
 41:       integer, parameter          :: image_cutoff = 5 39:       integer, parameter          :: image_cutoff = 5
 42:       LOGICAL          :: GTEST 40:       LOGICAL          :: GTEST
 43:  41: 
 44:       ! figure out how many lattice vectors to sum over 42:       !print *, 'hiya'
  43: 
  44:       !if (boxderivt) then
  45:       !   keep_angles = check_angles(box_params(4:6))
  46:       !   if (.not.keep_angles) then
  47:       !      print *, 'rejecting1'
  48:       !      call reject(energy, g)
  49:       !      return
  50:       !   endif
  51:       !endif
  52: 
  53:       !call build_H(H, H_grad, gtest)
  54:       !call inversematrix(H, H_inverse)
  55:       !call get_volume(vol)
  56:       !call get_reciplatvec(reciplatvec,reciplatvec_grad, .false.)
  57: 
 45:       newaldreal(:) = floor(ewaldrealc/box_params(1:3) + 0.5d0) 58:       newaldreal(:) = floor(ewaldrealc/box_params(1:3) + 0.5d0)
  59:       !newaldreal(1) = floor(ewaldrealc*dsqrt(sum(reciplatvec(1,:)**2))/(2.0d0*pi) + 0.5d0)
  60:       !newaldreal(2) = floor(ewaldrealc*dsqrt(sum(reciplatvec(2,:)**2))/(2.0d0*pi) + 0.5d0)
  61:       !newaldreal(3) = floor(ewaldrealc*dsqrt(sum(reciplatvec(3,:)**2))/(2.0d0*pi) + 0.5d0)
  62:       !print *, 'newaldreal: ', newaldreal(:3)
 46:  63: 
 47:       ! reject structure if would have to sum over more than five lattice vectors 
 48:       if (boxderivt) then 64:       if (boxderivt) then
 49:          if (.not. all(newaldreal.le.image_cutoff)) then 65:          if (.not. all(newaldreal.le.image_cutoff)) then
  66:             print *, 'rejecting2'
 50:             call reject(energy, g) 67:             call reject(energy, g)
 51:             return 68:             return
 52:          endif 69:          endif
 53:       endif 70:       endif
 54:  71: 
 55:       ! factorials 72:       ! factorials
 56:       FCT(1) = 1; FCT(2) = 2; FCT(3) = 6; FCT(4) = 24; FCT(5) = 120; FCT(6) = 720 73:       FCT(1) = 1; FCT(2) = 2; FCT(3) = 6; FCT(4) = 24; FCT(5) = 120; FCT(6) = 720
 57:       ! initialize energy values 74:       ! initialize energy values
 58:       ! energy1 is due to short-range anisotropic interactions 75:       ! energy1 is due to short-range anisotropic interactions
 59:       ! energy2 is due to damped dispersion 76:       ! energy2 is due to damped dispersion
 60:       ! energy3 is due to long-range electrostatics (computed using Ewald) 77:       ! energy3 is due to long-range electrostatics (computed using Ewald)
 61:       ENERGY = 0.D0; ENERGY1 = 0.D0; ENERGY2 = 0.D0; ENERGY3 = 0.D0 78:       ENERGY = 0.D0; ENERGY1 = 0.D0; ENERGY2 = 0.D0; ENERGY3 = 0.D0
 62:  79: 
 63:       ! initialize gradient if GTEST true 80:       ! initialize gradient if GTEST true
 64:       IF (GTEST) G(:) = 0.D0 81:       IF (GTEST) G(:) = 0.D0
 65:       IF (GTEST) G3C(:) = 0.D0 82:       IF (GTEST) G3C(:) = 0.D0
 66:  83: 
  84:       !print *, 'boxparams benzrigid: ', box_params(1:6)
  85:       !print *, 'coords benzrigid   : ', x(1:25)
  86: 
 67:       ! dj337: check if input coordinates are cartesian 87:       ! dj337: check if input coordinates are cartesian
 68:       ! assumes ATOMRIGIDCOORDT is correct 88:       ! assumes ATOMRIGIDCOORDT is correct
 69:       IF (ATOMRIGIDCOORDT) THEN ! if input is cartesian 89:       IF (ATOMRIGIDCOORDT) THEN ! if input is cartesian
  90:          !print *, 'converting...'
 70:          ! convert to rigidbody coordinates 91:          ! convert to rigidbody coordinates
 71:          XR(:) = 0.D0 92:          XR(:) = 0.D0
 72:          CALL TRANSFORMCTORIGID(X, XR) 93:          CALL TRANSFORMCTORIGID(X, XR)
 73:          if (boxderivt) then 94:          if (boxderivt) then
 74:             call frac2cart_rb_ortho(nrigidbody, xdum, xr) 95:             call frac2cart_rb_ortho(nrigidbody, xdum, xr)
  96:             !if (ortho) call frac2cart_rb_ortho(xdum, xr)
 75:             x(:) = xdum(:) 97:             x(:) = xdum(:)
 76:          else 98:          else
 77:             x(:) = xr(:) 99:             x(:) = xr(:)
 78:          endif100:          endif
 79:       ENDIF101:       ENDIF
 80: 102: 
 81:       EWALDREALC2 = EWALDREALC**2 ! real-space cutoff103:       !print *, 'coords benzrigid   : '
 104:       !do j1 = 1, 2*nrigidbody
 105:       !   print *, x(3*j1-2:3*j1)
 106:       !enddo
 107: 
 108:       !call build_H(H, H_grad, gtest)
 109:       !call inversematrix(H, H_inverse)
 110:       !call get_volume(vol)
 111:       !call get_reciplatvec(reciplatvec,reciplatvec_grad, .false.)
 112:       !if (boxderivt) then
 113:       !   ! compute v factor
 114:       !   c(:) = dcos(box_params(4:6))
 115:       !   s(:) = dsin(box_params(4:6))
 116:       !   v_fact = dsqrt(1.0d0 - c(1)**2-c(2)**2-c(3)**2 + 2.0d0*c(1)*c(2)*c(3))
 117:       !   dv_fact(1) = s(1)*(c(1) - c(2)*c(3))/v_fact
 118:       !   dv_fact(2) = s(2)*(c(2) - c(1)*c(3))/v_fact
 119:       !   dv_fact(3) = s(3)*(c(3) - c(1)*c(2))/v_fact
 120:       !endif
 121: 
 122:       EWALDREALC2 = EWALDREALC**2
 123: 
 124:       !print *, 'reciplatvec: ', reciplatvec(:3,:3)
 82: 125: 
 83:       ! OFFSET is number of CoM coords (3*NRIGIDBODY)126:       ! OFFSET is number of CoM coords (3*NRIGIDBODY)
 84:       OFFSET     = 3*NRIGIDBODY127:       OFFSET     = 3*NRIGIDBODY
 85: 128: 
 86:       ! Computing Cartesian coordinates for the system.  129:       ! Computing Cartesian coordinates for the system.  
 87:       DO J1 = 1, NRIGIDBODY130:       DO J1 = 1, NRIGIDBODY
 88: 131: 
 89:          J3 = 3*J1132:          J3 = 3*J1
 90:          J5 = OFFSET + J3133:          J5 = OFFSET + J3
 91:          ! CoM coords for rigid body J1134:          ! CoM coords for rigid body J1
119:                DE1(J4,:) = MATMUL(DRMI1(:,:),RBSTLA(J2,:))162:                DE1(J4,:) = MATMUL(DRMI1(:,:),RBSTLA(J2,:))
120:                DE2(J4,:) = MATMUL(DRMI2(:,:),RBSTLA(J2,:))163:                DE2(J4,:) = MATMUL(DRMI2(:,:),RBSTLA(J2,:))
121:                DE3(J4,:) = MATMUL(DRMI3(:,:),RBSTLA(J2,:))164:                DE3(J4,:) = MATMUL(DRMI3(:,:),RBSTLA(J2,:))
122: 165: 
123:             ENDIF166:             ENDIF
124: 167: 
125:          ENDDO168:          ENDDO
126: 169: 
127:       ENDDO170:       ENDDO
128: 171: 
 172:       !print *, 'cart coords benzrigid   : '
 173:       !do j1 = 1, natoms
 174:       !   print *, r(j1,:3)
 175:       !enddo
 176: 
129:       ! Now compute the actual potential.177:       ! Now compute the actual potential.
130:       ! loop over rigid bodies (A)178:       ! loop over rigid bodies (A)
131:       DO J1 = 1, NRIGIDBODY - 1179:       DO J1 = 1, NRIGIDBODY - 1
132: 180: 
133:          J3 = 3*J1181:          J3 = 3*J1
134:          J5 = OFFSET + J3182:          J5 = OFFSET + J3
135:          ! CoM coords for rigid body J1183:          ! CoM coords for rigid body J1
136:          RI(:)  = X(J3-2:J3)184:          RI(:)  = X(J3-2:J3)
137: 185: 
138:          ! loop over sites in the rigid body J1186:          ! loop over sites in the rigid body J1
154: 202: 
155:                   ! J8 is index for site J203:                   ! J8 is index for site J
156:                   J8     = MAXSITE*(J2-1) + J204:                   J8     = MAXSITE*(J2-1) + J
157:                   ! EJ is Z-axis for site J205:                   ! EJ is Z-axis for site J
158:                   EJ(:)  = E(J8,:)206:                   EJ(:)  = E(J8,:)
159:                   rr(:) = r(j7,:) - r(j8,:)207:                   rr(:) = r(j7,:) - r(j8,:)
160:                   ! minimum image convention208:                   ! minimum image convention
161:                   rssmin(1) = rr(1) - box_params(1)*anint(rr(1)/box_params(1))209:                   rssmin(1) = rr(1) - box_params(1)*anint(rr(1)/box_params(1))
162:                   rssmin(2) = rr(2) - box_params(2)*anint(rr(2)/box_params(2))210:                   rssmin(2) = rr(2) - box_params(2)*anint(rr(2)/box_params(2))
163:                   rssmin(3) = rr(3) - box_params(3)*anint(rr(3)/box_params(3))211:                   rssmin(3) = rr(3) - box_params(3)*anint(rr(3)/box_params(3))
 212:                   ! convert to fractional coordinates
 213:                   !rrfrac(:) = matmul(H_inverse, rr(:))
 214:                   !! minimum image convention
 215:                   !rssfracmin(1) = rrfrac(1) - anint(rrfrac(1))
 216:                   !rssfracmin(2) = rrfrac(2) - anint(rrfrac(2))
 217:                   !rssfracmin(3) = rrfrac(3) - anint(rrfrac(3))
164: 218: 
165:                   if (gtest.and.boxderivt) then219:                   if (gtest.and.boxderivt) then
166:                      ! get center of mass separation vector220:                      ! get center of mass separation vector
167:                      rrcom(:) = x(j3-2:j3) - x(j4-2:j4)221:                      rrcom(:) = x(j3-2:j3) - x(j4-2:j4)
168:                      ! minimum image convention222:                      ! minimum image convention
169:                      rrcommin(1) = rrcom(1) - box_params(1)*anint(rr(1)/box_params(1))223:                      rrcommin(1) = rrcom(1) - box_params(1)*anint(rr(1)/box_params(1))
170:                      rrcommin(2) = rrcom(2) - box_params(2)*anint(rr(2)/box_params(2))224:                      rrcommin(2) = rrcom(2) - box_params(2)*anint(rr(2)/box_params(2))
171:                      rrcommin(3) = rrcom(3) - box_params(3)*anint(rr(3)/box_params(3))225:                      rrcommin(3) = rrcom(3) - box_params(3)*anint(rr(3)/box_params(3))
 226:                      ! convert to fractional coordinates
 227:                      !rrcomfrac(:) = matmul(H_inverse, rrcom(:))
 228:                      !! minimum image convention
 229:                      !rcomfracmin(1) = rrcomfrac(1) - anint(rrfrac(1))
 230:                      !rcomfracmin(2) = rrcomfrac(2) - anint(rrfrac(2))
 231:                      !rcomfracmin(3) = rrcomfrac(3) - anint(rrfrac(3))
172:                   endif232:                   endif
173: 233: 
174:                   ! sum over lattice vectors 
175:                   do l = -newaldreal(1), newaldreal(1)234:                   do l = -newaldreal(1), newaldreal(1)
176:                   rss(1) = rssmin(1) + box_params(1)*l235:                   rss(1) = rssmin(1) + box_params(1)*l
177: 236: 
178:                      do m = -newaldreal(2), newaldreal(2)237:                      do m = -newaldreal(2), newaldreal(2)
179:                      rss(2) = rssmin(2) + box_params(2)*m238:                      rss(2) = rssmin(2) + box_params(2)*m
180: 239: 
181:                         do n = -newaldreal(3), newaldreal(3)240:                         do n = -newaldreal(3), newaldreal(3)
182:                         rss(3) = rssmin(3) + box_params(3)*n241:                         rss(3) = rssmin(3) + box_params(3)*n
183: 242: 
 243:                         ! convert to absolute coordinates
 244:                         !rss(:) = matmul(H, rssfrac(:))
 245: 
184:                         ! get COM vector246:                         ! get COM vector
185:                         if (gtest.and.boxderivt) then247:                         if (gtest.and.boxderivt) then
186:                            rcom(1) = rrcommin(1) + box_params(1)*l248:                            rcom(1) = rrcommin(1) + box_params(1)*l
187:                            rcom(2) = rrcommin(2) + box_params(2)*m249:                            rcom(2) = rrcommin(2) + box_params(2)*m
188:                            rcom(3) = rrcommin(3) + box_params(3)*n250:                            rcom(3) = rrcommin(3) + box_params(3)*n
189:                         endif251:                         endif
190:                      252:                      
191:                         R2     = DOT_PRODUCT(RSS(:),RSS(:))253:                         R2     = DOT_PRODUCT(RSS(:),RSS(:))
 254:                         !print *, 'r2: ', r2
192:                         ! check if distance within cutoff255:                         ! check if distance within cutoff
193:                         IF (R2 < EWALDREALC2) THEN256:                         IF (R2 < EWALDREALC2) THEN
 257:                            !print *, j7, j8
 258:                            !print *, 'r   : ', rss(:3)
 259:                            !print *, 'rmin: ', rssmin(1:3)
 260:                            !print *, 'rcom: ', rcom(1:3)
194:                            ! ABSRIJ is site-site separation between I and J261:                            ! ABSRIJ is site-site separation between I and J
195:                            ABSRIJ = DSQRT(R2)262:                            ABSRIJ = DSQRT(R2)
196:                            ! NR is unit site-site vector from sites I to J263:                            ! NR is unit site-site vector from sites I to J
197:                            NR(:)  = RSS(:)/ABSRIJ264:                            NR(:)  = RSS(:)/ABSRIJ
198:                            R2     = 1.D0/R2265:                            R2     = 1.D0/R2
199:                            R6     = R2*R2*R2266:                            R6     = R2*R2*R2
 267:                            !print *, 'absrij: ', absrij
200:          268:          
201:       !     CALCULATE THE DISPERSION DAMPING FACTOR269:       !     CALCULATE THE DISPERSION DAMPING FACTOR
202:          270:          
203:                            ! initialize sum for the damping function and vertical shift271:                            ! initialize sum for the damping function and vertical shift
204:                            DMPFCT = 1.D0272:                            DMPFCT = 1.D0
205:                            DMPFCT_SHIFT = 1.D0273:                            DMPFCT_SHIFT = 1.D0
206:                            ! initialize sum for the derivative of damping function274:                            ! initialize sum for the derivative of damping function
207:                            DDMPDR = B275:                            DDMPDR = B
208:          276:          
209:                            ! calculate sums277:                            ! calculate sums
279:                               RHOCC   = RHOCC0 + RHOCC10*(COSTA + COSTB) + RHOCC20*(1.5D0*COSTA*COSTA & 347:                               RHOCC   = RHOCC0 + RHOCC10*(COSTA + COSTB) + RHOCC20*(1.5D0*COSTA*COSTA & 
280:                                       + 1.5D0*COSTB*COSTB - 1.D0)348:                                       + 1.5D0*COSTB*COSTB - 1.D0)
281:                               ! ENERGY1 is energy due to short-range anisotropic interactions349:                               ! ENERGY1 is energy due to short-range anisotropic interactions
282:                               ! calculate vertical shift for first term350:                               ! calculate vertical shift for first term
283:                               EXPFCT  = KKJ*DEXP(-ALPHACC*(ABSRIJ - RHOCC))351:                               EXPFCT  = KKJ*DEXP(-ALPHACC*(ABSRIJ - RHOCC))
284:                               VSHIFT1 = KKJ*DEXP(-ALPHACC*(EWALDREALC - RHOCC))352:                               VSHIFT1 = KKJ*DEXP(-ALPHACC*(EWALDREALC - RHOCC))
285:                               ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1353:                               ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
286:                               ! ENERGY2 is energy due to damped dispersion354:                               ! ENERGY2 is energy due to damped dispersion
287:                               ! calculate vertical shift for second term355:                               ! calculate vertical shift for second term
288:                               VSHIFT2 = DC6CC*DMPFCT_SHIFT/(EWALDREALC**6)356:                               VSHIFT2 = DC6CC*DMPFCT_SHIFT/(EWALDREALC**6)
 357:                               !print *, 'energy: ', dc6cc*dmpfct*r6
289:                               ENERGY2 = ENERGY2 - DC6CC*DMPFCT*R6 + VSHIFT2358:                               ENERGY2 = ENERGY2 - DC6CC*DMPFCT*R6 + VSHIFT2
290:          359:          
291:                               IF (GTEST) THEN360:                               IF (GTEST) THEN
292:          361:          
293:                                  ! DVDR is derivative of dispersion damping factor energy with factor 1/Rab362:                                  ! DVDR is derivative of dispersion damping factor energy with factor 1/Rab
294:                                  DVDR    = 6.D0*DC6CC*R6*R2*DMPFCT - DC6CC*R6*DDMPDR 363:                                  DVDR    = 6.D0*DC6CC*R6*R2*DMPFCT - DC6CC*R6*DDMPDR 
 364:                                  !print *, 'grad: ', dvdr
295:                                  ! FRIJ is derivative of ENERGY1 wrt r_ij with factor 1/Rab365:                                  ! FRIJ is derivative of ENERGY1 wrt r_ij with factor 1/Rab
296:                                  FRIJ(:) = ALPHACC*EXPFCT*(-NR(:) + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADR(:) &366:                                  FRIJ(:) = ALPHACC*EXPFCT*(-NR(:) + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADR(:) &
297:                                          + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDR(:))367:                                          + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDR(:))
298:                                  ! TIJ is derivative of ENERGY1 wrt pi with factor 1/Rab368:                                  ! TIJ is derivative of ENERGY1 wrt pi with factor 1/Rab
299:                                  TIJ(:)  = ALPHACC*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADPI(:) &369:                                  TIJ(:)  = ALPHACC*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADPI(:) &
300:                                          + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDPI(:))370:                                          + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDPI(:))
301:                                  ! TJI is derivative of ENERGY1 wrt pj with factor 1/Rab371:                                  ! TJI is derivative of ENERGY1 wrt pj with factor 1/Rab
302:                                  TJI(:)  = ALPHACC*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADPJ(:) &372:                                  TJI(:)  = ALPHACC*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADPJ(:) &
303:                                          + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDPJ(:)) 373:                                          + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDPJ(:)) 
304:          374:          
306:          376:          
307:                            ! calculate if I and J are both hydorgens377:                            ! calculate if I and J are both hydorgens
308:                            ELSEIF (I > NCARBON .AND. J > NCARBON) THEN378:                            ELSEIF (I > NCARBON .AND. J > NCARBON) THEN
309:          379:          
310:                               RHOHH  = RHOHH0 + RHOHH10*(COSTA + COSTB) + RHOHH20*(1.5D0*COSTA*COSTA      &380:                               RHOHH  = RHOHH0 + RHOHH10*(COSTA + COSTB) + RHOHH20*(1.5D0*COSTA*COSTA      &
311:                                      + 1.5D0*COSTB*COSTB - 1.D0) 381:                                      + 1.5D0*COSTB*COSTB - 1.D0) 
312:                               EXPFCT  = KKJ*DEXP(-ALPHAHH*(ABSRIJ - RHOHH))382:                               EXPFCT  = KKJ*DEXP(-ALPHAHH*(ABSRIJ - RHOHH))
313:                               VSHIFT1 = KKJ*DEXP(-ALPHAHH*(EWALDREALC - RHOHH))383:                               VSHIFT1 = KKJ*DEXP(-ALPHAHH*(EWALDREALC - RHOHH))
314:                               ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1384:                               ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
315:                               VSHIFT2 = DC6HH*DMPFCT_SHIFT/(EWALDREALC**6)385:                               VSHIFT2 = DC6HH*DMPFCT_SHIFT/(EWALDREALC**6)
 386:                               !print *, 'energy: ', dc6hh*dmpfct*r6
316:                               ENERGY2 = ENERGY2 - DC6HH*DMPFCT*R6 + VSHIFT2387:                               ENERGY2 = ENERGY2 - DC6HH*DMPFCT*R6 + VSHIFT2
317:          388:          
318:                               IF (GTEST) THEN389:                               IF (GTEST) THEN
319:          390:          
320:                                  DVDR    = 6.D0*DC6HH*R6*R2*DMPFCT - DC6HH*R6*DDMPDR 391:                                  DVDR    = 6.D0*DC6HH*R6*R2*DMPFCT - DC6HH*R6*DDMPDR 
 392:                                  !print *, 'grad: ', dvdr
321:                                  FRIJ(:) = ALPHAHH*EXPFCT*(-NR(:) + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADR(:) &393:                                  FRIJ(:) = ALPHAHH*EXPFCT*(-NR(:) + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADR(:) &
322:                                          + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDR(:))394:                                          + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDR(:))
323:                                  TIJ(:)  = ALPHAHH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPI(:) &395:                                  TIJ(:)  = ALPHAHH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPI(:) &
324:                                          + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPI(:))396:                                          + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPI(:))
325:                                  TJI(:)  = ALPHAHH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPJ(:) &397:                                  TJI(:)  = ALPHAHH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPJ(:) &
326:                                          + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPJ(:))398:                                          + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPJ(:))
327:          399:          
328:                               ENDIF400:                               ENDIF
329:          401:          
330:                            ! calculate if I is carbon and J is hydrogen402:                            ! calculate if I is carbon and J is hydrogen
331:                            ELSE IF (I <= NCARBON .AND. J > NCARBON) THEN 403:                            ELSE IF (I <= NCARBON .AND. J > NCARBON) THEN 
332:          404:          
333:                               RHOCH  = RHOCH0 + RHOC10H*COSTA + RHOCH10*COSTB + RHOC20H*(1.5D0*COSTA*COSTA &405:                               RHOCH  = RHOCH0 + RHOC10H*COSTA + RHOCH10*COSTB + RHOC20H*(1.5D0*COSTA*COSTA &
334:                                      - 0.5D0) + RHOCH20*(1.5D0*COSTB*COSTB - 0.5D0)406:                                      - 0.5D0) + RHOCH20*(1.5D0*COSTB*COSTB - 0.5D0)
335:                               EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))407:                               EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))
336:                               VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))408:                               VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))
337:                               ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1409:                               ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
338:                               VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)410:                               VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)
 411:                               !print *, 'energy: ', dc6ch*dmpfct*r6
339:                               ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2412:                               ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2
340:          413:          
341:                               IF (GTEST) THEN414:                               IF (GTEST) THEN
342:                            415:                            
343:                                  DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 416:                                  DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 
 417:                                  !print *, 'grad: ', dvdr
344:                                  FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADR(:) &418:                                  FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADR(:) &
345:                                          + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDR(:))419:                                          + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDR(:))
346:                                  TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPI(:) &420:                                  TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPI(:) &
347:                                          + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPI(:))421:                                          + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPI(:))
348:                                  TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPJ(:) &422:                                  TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPJ(:) &
349:                                          + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPJ(:))423:                                          + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPJ(:))
350:          424:          
351:                               ENDIF425:                               ENDIF
352:          426:          
353:                            ELSE !IF(I > NCARBON .AND. J <= NCARBON) THEN427:                            ELSE !IF(I > NCARBON .AND. J <= NCARBON) THEN
354:          428:          
355:                               RHOCH  = RHOCH0 + RHOCH10*COSTA + RHOC10H*COSTB + RHOCH20*(1.5D0*COSTA*COSTA &429:                               RHOCH  = RHOCH0 + RHOCH10*COSTA + RHOC10H*COSTB + RHOCH20*(1.5D0*COSTA*COSTA &
356:                                      - 0.5D0) + RHOC20H*(1.5D0*COSTB*COSTB - 0.5D0)430:                                      - 0.5D0) + RHOC20H*(1.5D0*COSTB*COSTB - 0.5D0)
357:                               EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))431:                               EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))
358:                               VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))432:                               VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))
359:                               ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1433:                               ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
360:                               VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)434:                               VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)
 435:                               !print *, 'energy: ', dc6ch*dmpfct*r6
361:                               ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2436:                               ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2
362:          437:          
363:                               IF (GTEST) THEN438:                               IF (GTEST) THEN
364:          439:          
365:                                  DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 440:                                  DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 
 441:                                  !print *, 'grad: ', dvdr
366:                                  FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADR(:) &442:                                  FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADR(:) &
367:                                          + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDR(:))443:                                          + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDR(:))
368:                                  TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPI(:) &444:                                  TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPI(:) &
369:                                          + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPI(:))445:                                          + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPI(:))
370:                                  TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPJ(:) &446:                                  TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPJ(:) &
371:                                          + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPJ(:))447:                                          + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPJ(:))
372:          448:          
373:                               ENDIF449:                               ENDIF
374:          450:          
375:                            ENDIF451:                            ENDIF
376:          452:          
377:                            IF (GTEST) THEN453:                            IF (GTEST) THEN
378:          454:          
379:                               ! total gradient wrt CoM coords for rigid body J1455:                               ! total gradient wrt CoM coords for rigid body J1
380:                               G(J3-2:J3) = G(J3-2:J3) + DVDR*RSS(:) + FRIJ(:)456:                               G(J3-2:J3) = G(J3-2:J3) + DVDR*RSS(:) + FRIJ(:)
 457:                               !g(j3-2:j3) = g(j3-2:j3) + frij(:)
 458: 
 459:                               !box_paramsgrad(1:3) = box_paramsgrad(1:3) + (dvdr*rss(1:3)+frij(1:3))*rcom(1:3)/box_params(1:3)
381:                               ! total gradient wrt CoM coords for rigid body J2460:                               ! total gradient wrt CoM coords for rigid body J2
382:                               G(J4-2:J4) = G(J4-2:J4) - DVDR*RSS(:) - FRIJ(:)461:                               G(J4-2:J4) = G(J4-2:J4) - DVDR*RSS(:) - FRIJ(:)
 462:                               !g(j4-2:j4) = g(j4-2:j4) - frij(:)
383: 463: 
384:                               ! total gradient wrt AA coords for rigid body J1464:                               ! total gradient wrt AA coords for rigid body J1
385:                               G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPI(:) + TIJ(:)465:                               G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPI(:) + TIJ(:)
 466:                               !g(j5-2:j5) = g(j5-2:j5) + tij(:)
386:                               ! total gradient wrt AA coords for rigid body J2467:                               ! total gradient wrt AA coords for rigid body J2
387:                               G(J6-2:J6) = G(J6-2:J6) + DVDR*DRIJDPJ(:) + TJI(:)468:                               G(J6-2:J6) = G(J6-2:J6) + DVDR*DRIJDPJ(:) + TJI(:)
 469:                               !g(j6-2:j6) = g(j6-2:j6) + tji(:)
 470: 
 471:                               ! TODO: this is if orientation of rigid body depends on cell parameters
 472:                               !H_mat(:,:) = (-1.0d0/(vol*v_fact))*dv_fact(1)*H(:,:) + (1.0d0/vol)*H_grad(:,:,4)
 473:                               !rmatrix(:) = matmul(H_inverse, vol*(rss(:)-rcom(:)))
 474:                               !box_paramsgrad(4) = box_paramsgrad(4) + dvdr*dot_product(rss(1:3), matmul(H_mat, rmatrix(:)))
388: 475: 
389:                               ! gradients wrt cell lengths 
390:                               if (boxderivt) box_paramsgrad(1:3) = box_paramsgrad(1:3) + (dvdr*rss(1:3)+frij(1:3))*rcom(1:3)/box_params(1:3)476:                               if (boxderivt) box_paramsgrad(1:3) = box_paramsgrad(1:3) + (dvdr*rss(1:3)+frij(1:3))*rcom(1:3)/box_params(1:3)
 477:                               !box_paramsgrad(1) = box_paramsgrad(1) + dvdr*dot_product(rss(1:3), matmul(H_mat/vol, rmatrix(:)))
 478:                               !if (boxderivt) then
 479:                               !   do idx = 1, 6
 480:                               !      box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product((dvdr*rss(1:3) + frij(1:3)), matmul(H_grad(:,:,idx), rcomfrac(:)))
 481:                               !   enddo
 482:                               !endif
391: 483: 
392:                            ENDIF ! gtest484:                            ENDIF
393:       485:       
394:                         ENDIF ! within cutoff486:                         ENDIF
395: 487: 
396:                      enddo ! n488:                      enddo
397:                   enddo ! m489:                   enddo
398:                enddo ! l490:                enddo
399: 491: 
400:                ENDDO ! sites j492:                ENDDO
401: 493: 
402:             ENDDO ! rigid bodies J494:             ENDDO
403:  495:  
404:          ENDDO ! sites i496:          ENDDO
405: 497: 
406:       ENDDO ! rigid bodies I498:       ENDDO
407: 499: 
408: ! INCLUDE CONTRIBUTION OF RIGID BODY WITH PERIODIC IMAGE OF ITSELF500: ! INCLUDE CONTRIBUTION OF RIGID BODY WITH PERIODIC IMAGE OF ITSELF
409: 501: 
410:       ! loop over rigidbodies502:       ! loop over rigidbodies
411:       do j1 = 1, nrigidbody503:       do j1 = 1, nrigidbody
412:          j3 = 3*j1504:          j3 = 3*j1
413:          j5 = offset + j3505:          j5 = offset + j3
414:          ri(:) = x(j3-2:j3)506:          ri(:) = x(j3-2:j3)
415: 507: 
416:          ! loop over sites i508:          ! loop over sites i
417:          do i = 1, nsiteperbody(j1)509:          do i = 1, nsiteperbody(j1)
418:             j7 = maxsite*(j1-1) + i510:             j7 = maxsite*(j1-1) + i
419:             ei(:) = e(j7,:)511:             ei(:) = e(j7,:)
420: 512: 
421:             ! loop over sites j513:             ! loop over sites j
422:             do j = 1, nsiteperbody(j1)514:             do j = 1, nsiteperbody(j1)
423:                j8 = maxsite*(j1-1) + j515:                j8 = maxsite*(j1-1) + j
424:                ej(:) = e(j8,:)516:                ej(:) = e(j8,:)
425: 517: 
426:                ! site-site separation vector518:                !print *, 'new!!', j1, i, j
427:                rr(:) = r(j7,:) - r(j8,:)519:                rr(:) = r(j7,:) - r(j8,:)
 520:                !print *, j7, r(j7,:3)
 521:                !print *, j8, r(j8,:3)
 522:                !print *, 'rr: ', rr(:3)
 523:                ! convert to fractional
 524:                !rrfrac(:) = matmul(H_inverse, rr(:))
 525:                !print *, 'rrfrac: ', rrfrac(:3)
428: 526: 
429:                ! sum over lattice vectors527:                ! loop over boxes
430:                do l = -newaldreal(1), newaldreal(1)528:                do l = -newaldreal(1), newaldreal(1)
431:                   do m = -newaldreal(2), newaldreal(2)529:                   do m = -newaldreal(2), newaldreal(2)
432:                      do n = -newaldreal(3), newaldreal(3)530:                      do n = -newaldreal(3), newaldreal(3)
433: 531: 
434:                      ! if not in same rigid body 
435:                      if (.not.(l.eq.0.and.m.eq.0.and.n.eq.0)) then532:                      if (.not.(l.eq.0.and.m.eq.0.and.n.eq.0)) then
436: 533: 
437:                         rss(1) = rr(1) + box_params(1)*l534:                         rss(1) = rr(1) + box_params(1)*l
438:                         rss(2) = rr(2) + box_params(2)*m535:                         rss(2) = rr(2) + box_params(2)*m
439:                         rss(3) = rr(3) + box_params(3)*n536:                         rss(3) = rr(3) + box_params(3)*n
 537:                         !print *, l, m, n
 538:                         !print *, 'rssfrac: ', rssfrac(:3)
 539:                         !rss(:) = matmul(H, rssfrac(:))
440: 540: 
441:                         ! get COM vector 
442:                         if (gtest.and.boxderivt) then541:                         if (gtest.and.boxderivt) then
443:                            rcom(1) = box_params(1)*l542:                            rcom(1) = box_params(1)*l
444:                            rcom(2) = box_params(2)*m543:                            rcom(2) = box_params(2)*m
445:                            rcom(3) = box_params(3)*n544:                            rcom(3) = box_params(3)*n
446:                         endif545:                         endif
447: 546: 
448:                         r2 = dot_product(rss(:), rss(:))547:                         r2 = dot_product(rss(:), rss(:))
 548:                         !print *, 'rssmin2: ', rss(1:3)
449:                         if (r2 < ewaldrealc2) then549:                         if (r2 < ewaldrealc2) then
450: 550: 
451:                         ! absolute site-site distance551:                         !print *, 'rssmin3: ', rss(1:3)
 552:                         !print *, 'r2    : ', r2
 553:                         !print *, 'rcom  : ', rcom(1:3)
452:                         absrij = dsqrt(r2)554:                         absrij = dsqrt(r2)
 555:                         !print *, j7, j8
 556:                         !print *, absrij
453:                         nr(:) = rss(:)/absrij557:                         nr(:) = rss(:)/absrij
454:                         r2 = 1.d0/r2558:                         r2 = 1.d0/r2
455:                         r6 = r2*r2*r2559:                         r6 = r2*r2*r2
456: 560: 
457:                         ! CALCULATE DISPERSION DAMPING FACTOR561:                         ! CALCULATE DISPERSION DAMPING FACTOR
458: 562: 
459:                         ! initialize sum for the damping function and vertical shift563:                         ! initialize sum for the damping function and vertical shift
460:                         DMPFCT = 1.D0564:                         DMPFCT = 1.D0
461:                         DMPFCT_SHIFT = 1.D0565:                         DMPFCT_SHIFT = 1.D0
462:                         ! initialize sum for the derivative of damping function566:                         ! initialize sum for the derivative of damping function
534:                            RHOCC   = RHOCC0 + RHOCC10*(COSTA + COSTB) + RHOCC20*(1.5D0*COSTA*COSTA & 638:                            RHOCC   = RHOCC0 + RHOCC10*(COSTA + COSTB) + RHOCC20*(1.5D0*COSTA*COSTA & 
535:                                    + 1.5D0*COSTB*COSTB - 1.D0)639:                                    + 1.5D0*COSTB*COSTB - 1.D0)
536:                            ! ENERGY1 is energy due to short-range anisotropic interactions640:                            ! ENERGY1 is energy due to short-range anisotropic interactions
537:                            ! calculate vertical shift for first term641:                            ! calculate vertical shift for first term
538:                            EXPFCT  = KKJ*DEXP(-ALPHACC*(ABSRIJ - RHOCC))642:                            EXPFCT  = KKJ*DEXP(-ALPHACC*(ABSRIJ - RHOCC))
539:                            VSHIFT1 = KKJ*DEXP(-ALPHACC*(EWALDREALC - RHOCC))643:                            VSHIFT1 = KKJ*DEXP(-ALPHACC*(EWALDREALC - RHOCC))
540:                            ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1644:                            ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
541:                            ! ENERGY2 is energy due to damped dispersion645:                            ! ENERGY2 is energy due to damped dispersion
542:                            ! calculate vertical shift for second term646:                            ! calculate vertical shift for second term
543:                            VSHIFT2 = DC6CC*DMPFCT_SHIFT/(EWALDREALC**6)647:                            VSHIFT2 = DC6CC*DMPFCT_SHIFT/(EWALDREALC**6)
 648:                            !print *, 'energy: ', dc6cc*dmpfct*r6
544:                            ENERGY2 = ENERGY2 - DC6CC*DMPFCT*R6 + VSHIFT2649:                            ENERGY2 = ENERGY2 - DC6CC*DMPFCT*R6 + VSHIFT2
 650:                            !print *, 'vshift2     : ', vshift2
 651:                            !print *, 'contribution: ', -dc6cc*dmpfct*r6
 652:                            !print *, 'energy2     : ', energy2
545: 653: 
546:                            IF (GTEST) THEN654:                            IF (GTEST) THEN
547: 655: 
548:                               ! DVDR is derivative of dispersion damping factor energy with factor 1/Rab656:                               ! DVDR is derivative of dispersion damping factor energy with factor 1/Rab
549:                               DVDR    = 6.D0*DC6CC*R6*R2*DMPFCT - DC6CC*R6*DDMPDR 657:                               DVDR    = 6.D0*DC6CC*R6*R2*DMPFCT - DC6CC*R6*DDMPDR 
550:                               ! FRIJ is derivative of ENERGY1 wrt r_ij with factor 1/Rab658:                               !print *, 'dvdr        : ', dvdr
 659:                            !   !print *, 'grad: ', dvdr
 660:                            !   ! FRIJ is derivative of ENERGY1 wrt r_ij with factor 1/Rab
551:                               FRIJ(:) = ALPHACC*EXPFCT*(-NR(:) + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADR(:) &661:                               FRIJ(:) = ALPHACC*EXPFCT*(-NR(:) + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADR(:) &
552:                                       + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDR(:))662:                                       + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDR(:))
553:                               ! TIJ is derivative of ENERGY1 wrt pi with factor 1/Rab663:                               ! TIJ is derivative of ENERGY1 wrt pi with factor 1/Rab
554:                               TIJ(:)  = ALPHACC*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADPI(:) &664:                               TIJ(:)  = ALPHACC*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADPI(:) &
555:                                       + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDPI(:))665:                                       + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDPI(:))
556:                               ! TJI is derivative of ENERGY1 wrt pj with factor 1/Rab666:                            !   ! TJI is derivative of ENERGY1 wrt pj with factor 1/Rab
557:                               TJI(:)  = ALPHACC*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADPJ(:) &667:                               TJI(:)  = ALPHACC*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADPJ(:) &
558:                                       + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDPJ(:)) 668:                                       + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDPJ(:)) 
559: 669: 
560:                            ENDIF670:                            ENDIF
561: 671: 
562:                         ! calculate if I and J are both hydorgens672:                         ! calculate if I and J are both hydorgens
563:                         ELSEIF (I > NCARBON .AND. J > NCARBON) THEN673:                         ELSEIF (I > NCARBON .AND. J > NCARBON) THEN
564: 674: 
565:                            RHOHH  = RHOHH0 + RHOHH10*(COSTA + COSTB) + RHOHH20*(1.5D0*COSTA*COSTA      &675:                            RHOHH  = RHOHH0 + RHOHH10*(COSTA + COSTB) + RHOHH20*(1.5D0*COSTA*COSTA      &
566:                                   + 1.5D0*COSTB*COSTB - 1.D0)676:                                   + 1.5D0*COSTB*COSTB - 1.D0)
567:                            EXPFCT  = KKJ*DEXP(-ALPHAHH*(ABSRIJ - RHOHH))677:                            EXPFCT  = KKJ*DEXP(-ALPHAHH*(ABSRIJ - RHOHH))
568:                            VSHIFT1 = KKJ*DEXP(-ALPHAHH*(EWALDREALC - RHOHH))678:                            VSHIFT1 = KKJ*DEXP(-ALPHAHH*(EWALDREALC - RHOHH))
569:                            ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1679:                            ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
570:                            VSHIFT2 = DC6HH*DMPFCT_SHIFT/(EWALDREALC**6)680:                            VSHIFT2 = DC6HH*DMPFCT_SHIFT/(EWALDREALC**6)
 681:                            !print *, 'energy: ', dc6hh*dmpfct*r6
571:                            ENERGY2 = ENERGY2 - DC6HH*DMPFCT*R6 + VSHIFT2682:                            ENERGY2 = ENERGY2 - DC6HH*DMPFCT*R6 + VSHIFT2
572: 683: 
573:                            IF (GTEST) THEN684:                            IF (GTEST) THEN
574: 685: 
575:                               DVDR    = 6.D0*DC6HH*R6*R2*DMPFCT - DC6HH*R6*DDMPDR 686:                               DVDR    = 6.D0*DC6HH*R6*R2*DMPFCT - DC6HH*R6*DDMPDR 
 687:                            !   !print *, 'grad: ', dvdr
576:                               FRIJ(:) = ALPHAHH*EXPFCT*(-NR(:) + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADR(:) &688:                               FRIJ(:) = ALPHAHH*EXPFCT*(-NR(:) + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADR(:) &
577:                                       + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDR(:))689:                                       + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDR(:))
578:                               TIJ(:)  = ALPHAHH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPI(:) &690:                               TIJ(:)  = ALPHAHH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPI(:) &
579:                                       + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPI(:))691:                                       + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPI(:))
580:                               TJI(:)  = ALPHAHH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPJ(:) &692:                               TJI(:)  = ALPHAHH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPJ(:) &
581:                                       + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPJ(:))693:                                       + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPJ(:))
582: 694: 
583:                            ENDIF695:                            ENDIF
584: 696: 
585:                         ! calculate if I is carbon and J is hydrogen697:                         ! calculate if I is carbon and J is hydrogen
586:                         ELSE IF (I <= NCARBON .AND. J > NCARBON) THEN 698:                         ELSE IF (I <= NCARBON .AND. J > NCARBON) THEN 
587: 699: 
588:                            RHOCH  = RHOCH0 + RHOC10H*COSTA + RHOCH10*COSTB + RHOC20H*(1.5D0*COSTA*COSTA &700:                            RHOCH  = RHOCH0 + RHOC10H*COSTA + RHOCH10*COSTB + RHOC20H*(1.5D0*COSTA*COSTA &
589:                                   - 0.5D0) + RHOCH20*(1.5D0*COSTB*COSTB - 0.5D0)701:                                   - 0.5D0) + RHOCH20*(1.5D0*COSTB*COSTB - 0.5D0)
590:                            EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))702:                            EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))
591:                            VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))703:                            VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))
592:                            ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1704:                            ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
593:                            VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)705:                            VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)
 706:                            !print *, 'energy: ', dc6ch*dmpfct*r6
594:                            ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2707:                            ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2
595: 708: 
596:                            IF (GTEST) THEN709:                            IF (GTEST) THEN
597: 710: 
598:                               DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 711:                               DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 
 712:                            !   !print *, 'grad: ', dvdr
599:                               FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADR(:) &713:                               FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADR(:) &
600:                                       + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDR(:))714:                                       + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDR(:))
601:                               TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPI(:) &715:                               TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPI(:) &
602:                                       + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPI(:))716:                                       + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPI(:))
603:                               TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPJ(:) &717:                               TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPJ(:) &
604:                                       + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPJ(:))718:                                       + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPJ(:))
605: 719: 
606:                            ENDIF720:                            ENDIF
607: 721: 
608:                         ELSE !IF(I > NCARBON .AND. J <= NCARBON) THEN722:                         ELSE !IF(I > NCARBON .AND. J <= NCARBON) THEN
609: 723: 
610:                            RHOCH  = RHOCH0 + RHOCH10*COSTA + RHOC10H*COSTB + RHOCH20*(1.5D0*COSTA*COSTA &724:                            RHOCH  = RHOCH0 + RHOCH10*COSTA + RHOC10H*COSTB + RHOCH20*(1.5D0*COSTA*COSTA &
611:                                   - 0.5D0) + RHOC20H*(1.5D0*COSTB*COSTB - 0.5D0)725:                                   - 0.5D0) + RHOC20H*(1.5D0*COSTB*COSTB - 0.5D0)
612:                            EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))726:                            EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))
613:                            VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))727:                            VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))
614:                            ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1728:                            ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
615:                            VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)729:                            VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)
 730:                            !print *, 'energy: ', dc6ch*dmpfct*r6
616:                            ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2731:                            ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2
617: 732: 
618:                            IF (GTEST) THEN733:                            IF (GTEST) THEN
619: 734: 
620:                               DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 735:                               DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 
 736:                            !   !print *, 'grad: ', dvdr
621:                               FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADR(:) &737:                               FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADR(:) &
622:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDR(:))738:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDR(:))
623:                               TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPI(:) &739:                               TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPI(:) &
624:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPI(:))740:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPI(:))
625:                               TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPJ(:) &741:                               TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPJ(:) &
626:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPJ(:))742:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPJ(:))
627: 743: 
628:                            ENDIF744:                            ENDIF
629: 745: 
630:                         ENDIF746:                         ENDIF
631: 747: 
632: 748: 
633:                         IF (GTEST) THEN749:                         IF (GTEST) THEN
634: 750: 
635:                            ! total gradient wrt AA coords for rigid body J1751:                            ! total gradient wrt AA coords for rigid body J1
636:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPI(:) + TIJ(:)752:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPI(:) + TIJ(:)
 753:                            !g(j5-2:j5) = g(j5-2:j5) + tij(:)
637:                            ! total gradient wrt AA coords for rigid body J2754:                            ! total gradient wrt AA coords for rigid body J2
638:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPJ(:) + TJI(:)755:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPJ(:) + TJI(:)
 756:                            !g(j5-2:j5) = g(j5-2:j5) + tji(:)
 757:                            !print *, 'dispersion:'
 758:                            !print *, dvdr*drijdpi(:3)
 759:                            !print *, dvdr*drijdpj(:3)
 760:                            !print *, 'energy2: ', energy2
 761:                            !print *, 'anisotropic:'
 762:                            !print *, tij(:3)
 763:                            !print *, tji(:3)
639: 764: 
640:                            ! gradient wrt cell lengths 
641:                            if (boxderivt) box_paramsgrad(1:3) = box_paramsgrad(1:3) + (dvdr*rss(1:3)+frij(1:3))*rcom(1:3)/box_params(1:3)765:                            if (boxderivt) box_paramsgrad(1:3) = box_paramsgrad(1:3) + (dvdr*rss(1:3)+frij(1:3))*rcom(1:3)/box_params(1:3)
 766:                            !if (boxderivt) then
 767:                            !   do idx = 1, 6
 768:                            !         box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product((dvdr*rss(1:3) + frij(1:3)), matmul(H_grad(:,:,idx), rcomfrac(:)))
 769:                            !         !box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product(dvdr*rss(1:3), matmul(H_grad(:,:,idx), rcomfrac(:)))
 770:                            !   enddo
 771:                            !endif
642: 772: 
643:                         ENDIF ! gtest773:                         ENDIF
644:                         endif ! central box774:                         endif
645:                     endif ! within cutoff775:                     endif
646:                   enddo ! n776:                   enddo
647:                enddo ! m777:                enddo
648:             enddo ! l778:             enddo
649:             enddo ! sites j779:             enddo
650:          enddo ! sites i780:          enddo
651:       enddo ! rigid bodies781:       enddo
652: 782: 
653:       ! convert to cartesian coordinates783:       ! convert to cartesian coordinates
654:       XC(:) = 0.D0784:       XC(:) = 0.D0
655:       if (boxderivt) then785:       if (boxderivt) then
656:          xdum(:) = x(:)786:          xdum(:) = x(:)
657:          call cart2frac_rb_ortho(nrigidbody, xdum, x)787:          call cart2frac_rb_ortho(nrigidbody, xdum, x)
658:       endif788:       endif
659:       CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XC, X)789:       CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XC, X)
660:       ! restore cartesian rigid body coordinates790:       ! restore cartesian rigid body coordinates
661:       if (boxderivt) x(:) = xdum(:)791:       if (boxderivt) x(:) = xdum(:)
662: 792: 
663:       ! ENERGY3 and G3 are energy and gradient due to electrostatics793: !      !!! ENERGY3 and G3 are energy and gradient due to electrostatics
664:       ! computed using Ewald summation794: !      !!! computed using Ewald summation
665:       CALL EWALDSUM(1, XC, G3C, ENERGY3, GTEST)795:       CALL EWALDSUM(1, XC, G3C, ENERGY3, GTEST)
666: 796: !
667:       ! convert Ewald contribution of gradient to rigidbody coordinates797: !      !!! convert Ewald contribution of gradient to rigidbody coordinates
668:       IF (GTEST) G3(:) = 0.D0798:       IF (GTEST) G3(:) = 0.D0
669:       CALL TRANSFORMGRAD(G3C, X, G3)799:       CALL TRANSFORMGRAD(G3C, X, G3)
 800: !!      !print *, 'energy2: ', energy2
 801: 
 802:       !energy = (energy1+energy2)*2625.499d0
 803:       !if (gtest) g(:) = g(:)*2625.499d0
 804:       !energy = energy2*2625.499d0
 805:       !energy = (energy1+energy3)*2625.499d0
 806:       !if (gtest) g(:) = g3(:)*2625.499d0
 807:       !ENERGY = (ENERGY1 + ENERGY2 + ENERGY3)*2625.499D0 
 808:       !IF (GTEST) G(:) = (G(:) + G3(:))*2625.499D0
 809:       !if (gtest) box_paramsgrad(1:6) = box_paramsgrad(1:6)*2625.499D0
 810:       !print *, 'box_paramsgrad: ', box_paramsgrad(1:3)
 811:       !print *, 'energy: ', energy2
670: 812: 
671:       ! dj337: if input was cartesian, convert back to cartesian813:       ! dj337: if input was cartesian, convert back to cartesian
672:       ! assumes ATOMRIGIDCOORDT is correct814:       ! assumes ATOMRIGIDCOORDT is correct
673:       IF (ATOMRIGIDCOORDT) THEN815:       IF (ATOMRIGIDCOORDT) THEN
674: 816: 
675:          ! convert to cartesian coordinates817:          ! convert to cartesian coordinates
676:          if (boxderivt) then818:          if (boxderivt) then
 819:             !print *, 'coords 708: ', x(1:3*natoms)
677:             xdum(:) = x(:)820:             xdum(:) = x(:)
678:             call cart2frac_rb_ortho(nrigidbody, xdum, x)821:             call cart2frac_rb_ortho(nrigidbody, xdum, x)
 822:             !if (ortho) call cart2frac_rb_ortho(xdum, x)
 823:             !print *, 'coords 711: ', x(1:3*natoms)
679:          endif824:          endif
680:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XR, X)825:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XR, X)
681:          X(:) = XR(:)826:          X(:) = XR(:)
682:       ENDIF827:       ENDIF
683: 828: 
684:       ! sum energies / gradients and convert to kJ/mol829:       !if (boxderivt) call constrain_volume(v_fact, dv_fact, energy1, box_paramsgrad(4:6), gtest)
 830: 
685:       ENERGY = (ENERGY1 + ENERGY2 + ENERGY3)*2625.499D0831:       ENERGY = (ENERGY1 + ENERGY2 + ENERGY3)*2625.499D0
686:       IF (GTEST) G(:) = (G(:) + G3(:))*2625.499D0832:       IF (GTEST) G(:) = (G(:) + G3(:))*2625.499D0
687:       if (gtest) box_paramsgrad(1:3) = box_paramsgrad(1:3)*2625.499D0833:       if (gtest) box_paramsgrad(1:3) = box_paramsgrad(1:3)*2625.499D0
 834:       !print *, 'coords benzrigid    : ', x(1:3*natoms)
 835:       !print *, 'box params benzrigid: ', box_params(1:6)
 836: 
 837:       !print *, 'energy: ', energy!1, energy2, energy3
 838:       !print *, 'grad: ', g(1)
 839:       !!print *, 'grad  : '
 840:       !do j1 = 1, 2*nrigidbody
 841:       !   print *, x(3*j1-2:3*j1)
 842:       !enddo
 843: 
 844:       !print *, 'contri: ', energy1, energy2, energy3
 845:       !print *, 'energy: ', energy !, box_params(1:6)
 846:       !print *, 'g1    : ', g(1:3)
 847:       !stop
688: 848: 
689:       END SUBROUTINE BENZGENRIGIDEWALD_ORTHO849:       END SUBROUTINE BENZGENRIGIDEWALD_ORTHO
690: 850: 
691: !     ----------------------------------------------------------------------------------------------851: !     ----------------------------------------------------------------------------------------------


r33132/box_derivatives.f90 2017-08-07 15:31:02.781686129 +0100 r33131/box_derivatives.f90 2017-08-07 15:31:05.713725136 +0100
  1: module box_derivatives  1: module box_derivatives
  2: use commons, only: natoms, box_params, box_paramsgrad  2: use commons, only: natoms, box_params, box_paramsgrad
  3: use genrigid, only: degfreedoms, nrigidbody  3: use genrigid, only: degfreedoms, nrigidbody
  4: use cartdist  4: use cartdist
  5:   5: 
  6: implicit none  6: implicit none
  7:   7: 
  8: public :: check_angles  8: public :: check_angles
  9:   9: 
  10: ! TODO
 10: ! dj337: module to convert gradient from absolute to  11: ! dj337: module to convert gradient from absolute to 
 11: ! fractional coordinates to perform basin-hopping on 12: ! fractional coordinates to perform basin-hopping on
 12: ! periodic systems with varying lattice parameters 13: ! periodic systems with varying lattice parameters
 13:  14: 
 14: contains 15: contains
 15:  16: 
 16: ! ----------------------------------------------------------------------------------- 17: ! -----------------------------------------------------------------------------------
 17: ! VARIABLES 18: ! VARIABLES
 18: ! x: atomic positions in absolute coordinates 19: ! x: atomic positions in absolute coordinates
 19: ! xfrac: atomic positions in fractional coordinates 20: ! xfrac: atomic positions in fractional coordinates
256:           ! convert gradient wrt COM positions from absolute to fractional257:           ! convert gradient wrt COM positions from absolute to fractional
257:           gradrfrac(j3-2:j3) = matmul(gradr(j3-2:j3), H)258:           gradrfrac(j3-2:j3) = matmul(gradr(j3-2:j3), H)
258:        enddo259:        enddo
259: 260: 
260:        ! gradient wrt AA coordinates are unchanged261:        ! gradient wrt AA coordinates are unchanged
261:        gradrfrac(3*nrigidbody+1:degfreedoms) = gradr(3*nrigidbody+1:degfreedoms)262:        gradrfrac(3*nrigidbody+1:degfreedoms) = gradr(3*nrigidbody+1:degfreedoms)
262: 263: 
263:        end subroutine boxderiv_rb_tri264:        end subroutine boxderiv_rb_tri
264: 265: 
265: ! -----------------------------------------------------------------------------------266: ! -----------------------------------------------------------------------------------
266: ! TAKES A BASIN-HOPPING STEP by making uniformly random changes to the cell lengths267: ! TODO:
267: ! and angles 
268:  
269:        subroutine bd_takestep(np)268:        subroutine bd_takestep(np)
270: 269: 
271:        use commons, only: ortho, box_params270:        use commons, only: ortho, box_params
272:        use vec3, only: vec_random271:        use vec3, only: vec_random
273: 272: 
274:        integer, intent(in)         :: np273:        integer, intent(in)         :: np
275:        double precision            :: new_angles(3)274:        double precision            :: new_angles(3)
276:        double precision, parameter :: max_length_step = 0.3d0275:        double precision, parameter :: max_length_step = 0.3d0
277:        double precision, parameter :: max_angle_step = 0.1d0276:        double precision, parameter :: max_angle_step = 0.1d0
278: 277: 
 278:        !print *, 'taking step!'
 279: 
279:        ! generate new box lengths280:        ! generate new box lengths
280:        box_params(1:3) = box_params(1:3) + vec_random() * max_length_step281:        box_params(1:3) = box_params(1:3) + vec_random() * max_length_step
281:        ! if triclinic282:        ! if triclinic
282:        if (.not.(ortho)) then283:        if (.not.(ortho)) then
283:           new_angles(:) = box_params(4:6) + vec_random() * max_angle_step284:           new_angles(:) = box_params(4:6) + vec_random() * max_angle_step
284:           ! check to make sure combination of new angles is valid 
285:           do while (.not.check_angles(new_angles(:)))285:           do while (.not.check_angles(new_angles(:)))
 286:              !print *, 'bad angles'
286:              new_angles(:) = box_params(4:6) + vec_random() * max_angle_step287:              new_angles(:) = box_params(4:6) + vec_random() * max_angle_step
287:           enddo288:           enddo
288:           box_params(4:6) = new_angles(:)289:           box_params(4:6) = new_angles(:)
 290:              ! generate new box angles
 291:              !new_angles(:) = box_params(4:6) + vec_random() * max_angle_step
 292:              ! see whether new angles are valid
 293:              !if (check_angles(new_angles(:))) then
 294:                 
 295:                 !box_params(4:6) = new_angles(:)
 296:                 !exit
 297:              !endif ! check angles
 298:           !enddo
289:        endif ! triclinic299:        endif ! triclinic
290: 300: 
291:        end subroutine bd_takestep301:        end subroutine bd_takestep
292: 302: 
293: ! -----------------------------------------------------------------------------------303: ! -----------------------------------------------------------------------------------
294: ! CHECKS SET OF TRICLINIC CELL ANGLES to make sure they are valid. Non-valid set of angles304: ! TODO:
295: ! corresponds to a cell with zero or imaginary volume. Function returns True if the set 
296: ! of angles meets the criteria for valid cell angles. 
297: 305: 
298:        pure logical function check_angles(angles)306:        pure logical function check_angles(angles)
299: 307: 
300:        implicit none308:        implicit none
301: 309: 
302:        double precision, intent(in) :: angles(3)310:        double precision, intent(in) :: angles(3)
303:        double precision             :: sums(4)311:        double precision             :: sums(4)
304:        double precision, parameter   :: pi = 3.141592654d0312:        double precision, parameter   :: pi = 3.141592654d0
305: 313: 
306:        ! calculate necessary sums314:        ! calculate necessary sums
307:        sums(1) =  angles(1) + angles(2) + angles(3)315:        sums(1) =  angles(1) + angles(2) + angles(3)
308:        sums(2) = -angles(1) + angles(2) + angles(3)316:        sums(2) = -angles(1) + angles(2) + angles(3)
309:        sums(3) =  angles(1) - angles(2) + angles(3)317:        sums(3) =  angles(1) - angles(2) + angles(3)
310:        sums(4) =  angles(1) + angles(2) - angles(3)318:        sums(4) =  angles(1) + angles(2) - angles(3)
311: 319: 
312:        ! check all sums are between 0 and 2pi and all angles between 0 and pi320:        ! check all sums are between 0 and 2pi and all angles between 0 and pi
313:        check_angles = all(sums.gt.0.0d0).and.all(sums.lt.2*pi).and.all(angles.gt.0.0d0).and.all(angles.lt.pi)321:        check_angles = all(sums.gt.0.0d0).and.all(sums.lt.2*pi).and.all(angles.gt.0.0d0).and.all(angles.lt.pi)
314:        end function check_angles322:        end function check_angles
315: 323: 
316: ! -----------------------------------------------------------------------------------324: ! -----------------------------------------------------------------------------------
317: ! REJECTS structures that are invalid. Tricks the LBFGS minimizer by returning a very 
318: ! large energy and a very small gradient, so will the quench will immediately be  
319: ! considered converged but the structure will never be saved as a low-energy minimum. 
320: 325: 
321:        subroutine reject(energy, grad)326:        subroutine reject(energy, grad)
322: 327: 
323:        implicit none328:        implicit none
324: 329: 
325:        double precision, intent(out) :: energy, grad(3*natoms)330:        double precision, intent(out) :: energy, grad(3*natoms)
326: 331: 
327:        energy = 1.0d20332:        energy = 1.0d20
328:        grad(:) = 1.0d-20333:        grad(:) = 1.0d-20
329: 334: 
330:        end subroutine reject335:        end subroutine reject
331: 336: 
332: ! -----------------------------------------------------------------------------------337: ! -----------------------------------------------------------------------------------
333: ! ADDS WCA-STYLE REPULSION term to the energy and gradients to repel the cell away 
334: ! from having zero volume. 
335: 338: 
336:        subroutine constrain_volume(v, v_deriv, energy, grad_angles, gtest)339:        subroutine constrain_volume(v, v_deriv, energy, grad_angles, gtest)
337: 340: 
338:        implicit none341:        implicit none
339: 342: 
340:        double precision, intent(in)    :: v, v_deriv(3)343:        double precision, intent(in)    :: v, v_deriv(3)
341:        double precision, intent(inout) :: energy, grad_angles(3)344:        double precision, intent(inout) :: energy, grad_angles(3)
342:        logical, intent(in)             :: gtest345:        logical, intent(in)             :: gtest
343:        double precision, parameter     :: v_sigma = 3.0d-1346:        double precision, parameter     :: v_sigma = 3.0d-1
344:        double precision, parameter     :: v_eps = 1.0d-3347:        double precision, parameter     :: v_eps = 1.0d-3
348:           energy = energy + 4.0d0*v_eps*((v_sigma/v)**12 - (v_sigma/v)**6) + v_eps351:           energy = energy + 4.0d0*v_eps*((v_sigma/v)**12 - (v_sigma/v)**6) + v_eps
349: 352: 
350:           if (gtest) then353:           if (gtest) then
351:              ! add gradient contribution354:              ! add gradient contribution
352:              grad_angles(:) = grad_angles(:) + 24.0d0*v_eps/v_sigma*((v_sigma/v)**7 - 2.0d0*(v_sigma/v)**13)*v_deriv(:)355:              grad_angles(:) = grad_angles(:) + 24.0d0*v_eps/v_sigma*((v_sigma/v)**7 - 2.0d0*(v_sigma/v)**13)*v_deriv(:)
353:           endif356:           endif
354:        endif357:        endif
355: 358: 
356:        end subroutine constrain_volume359:        end subroutine constrain_volume
357: 360: 
358: ! ---------------------------------------------------------------------------------361: ! TODO: REJECT subroutine
 362: 
359: ! 363: ! 
360: ! Rotates all rigid bodies after a step is taken in the cell parameters.364: ! Rotates all rigid bodies after a step is taken in the cell parameters.
361: ! VARIABLES365: ! VARIABLES
362: ! box_paramsold: unit cell lengths and angles from before the step was taken366: ! box_paramsold: unit cell lengths and angles from before the step was taken
363: ! TODO: not sure this is working properly (or that the equations are even valid...)367: ! TODO: not sure this is working properly (or that the equations are even valid...)
364: 368: 
365: 369: 
366: !       subroutine rotate_bodies(box_paramsold, xrfrac)370: !       subroutine rotate_bodies(box_paramsold, xrfrac)
367: !371: !
368: !       use genrigid, only: transformctorigid, sitesrigidbody, maxsite, nsiteperbody, inversematrix372: !       use genrigid, only: transformctorigid, sitesrigidbody, maxsite, nsiteperbody, inversematrix


r33132/cartdist.f90 2017-08-07 15:31:03.001689060 +0100 r33131/cartdist.f90 2017-08-07 15:31:05.933728061 +0100
  1: module cartdist  1: module cartdist
  2: use commons, only: natoms, box_params, box_paramsgrad  2: use commons, only: natoms, box_params, box_paramsgrad
  3:   3: 
  4: implicit none  4: implicit none
  5:   5: 
  6: ! dj337: util modules for converting between absolute and fractional coordinate  6: ! TODO
  7: ! systems, building the H matrix, and computing the reciprocal lattice vectors  7: ! dj337: util modules
  8:   8: 
  9: contains  9: contains
 10:  10: 
 11: ! ----------------------------------------------------------------------------------- 11: ! -----------------------------------------------------------------------------------
 12: ! VARIABLES 12: ! VARIABLES
 13: ! x: atomic positions in absolute coordinates 13: ! x: atomic positions in absolute coordinates
 14: ! xfrac: atomic positions in fractional coordinates 14: ! xfrac: atomic positions in fractional coordinates
 15:  15: 
 16: ! xr: absolute rigid body (COM+AA) coordinates 16: ! xr: absolute rigid body (COM+AA) coordinates
 17: ! xrfrac: fractional rigid body coordinates 17: ! xrfrac: fractional rigid body coordinates
236:       logical, intent(in)                    :: gtest236:       logical, intent(in)                    :: gtest
237:       !double precision, intent(in), optional :: box_parameters237:       !double precision, intent(in), optional :: box_parameters
238:       double precision                       :: box_lengths(3), box_angles(3)238:       double precision                       :: box_lengths(3), box_angles(3)
239:       double precision                       :: c(3), s(3), v239:       double precision                       :: c(3), s(3), v
240: 240: 
241:       H(:,:) = 0.0d0241:       H(:,:) = 0.0d0
242:       H_grad(:,:,:) = 0.0d0242:       H_grad(:,:,:) = 0.0d0
243:       box_lengths(:) = box_params(1:3)243:       box_lengths(:) = box_params(1:3)
244:       box_angles(:) = box_params(4:6)244:       box_angles(:) = box_params(4:6)
245: 245: 
 246:       !if (present(box_parameters)) then
 247:       !   box_lengths(:) = box_parameters(1:3)
 248:       !   box_angles(:) = box_parameters(4:6)
 249:       !else
 250:       !   box_lengths(:) = box_params(1:3)
 251:       !   box_angles(:) = box_params(4:6)
 252:       !endif
 253: 
246:       ! cosine of the angles254:       ! cosine of the angles
247:       c(:) = dcos(box_angles(:))255:       c(:) = dcos(box_angles(:))
248:       ! sine of the angles256:       ! sine of the angles
249:       s(:) = dsin(box_angles(:))257:       s(:) = dsin(box_angles(:))
250:       ! factor that is related to the volume (but not quite volume)258:       ! factor that is related to the volume (but not quite volume)
251:       v = dsqrt(1.0d0 - c(1)**2 - c(2)**2 - c(3)**2 + 2.0d0*c(1)*c(2)*c(3))259:       v = dsqrt(1.0d0 - c(1)**2 - c(2)**2 - c(3)**2 + 2.0d0*c(1)*c(2)*c(3))
252: 260: 
253:       ! define the H transformation matrix261:       ! define the H transformation matrix
254:       ! first row of matrix262:       ! first row of matrix
255:       H(1,1) = box_lengths(1)263:       H(1,1) = box_lengths(1)
301: 309: 
302:        implicit none310:        implicit none
303: 311: 
304:        !double precision, intent(in), optional :: box_parameters(6)312:        !double precision, intent(in), optional :: box_parameters(6)
305:        double precision, intent(out)          :: vol313:        double precision, intent(out)          :: vol
306:        double precision                       :: box_lengths(3), box_angles(3), c(3)314:        double precision                       :: box_lengths(3), box_angles(3), c(3)
307: 315: 
308:        box_lengths(:) = box_params(1:3)316:        box_lengths(:) = box_params(1:3)
309:        box_angles(:) = box_params(4:6)317:        box_angles(:) = box_params(4:6)
310: 318: 
 319:        !if (present(box_parameters)) then
 320:        !   box_lengths(:) = box_parameters(1:3)
 321:        !   box_angles(:) = box_parameters(4:6)
 322:        !else
 323:        !   box_lengths(:) = box_params(1:3)
 324:        !   box_angles(:) = box_params(4:6)
 325:        !endif
 326: 
311:        vol = box_lengths(1)*box_lengths(2)*box_lengths(3)327:        vol = box_lengths(1)*box_lengths(2)*box_lengths(3)
312:        if (.not.ortho) then328:        if (.not.ortho) then
313:           c(:) = dcos(box_angles(:))329:           c(:) = dcos(box_angles(:))
314:           vol = vol * dsqrt(1.0d0 - c(1)**2 - c(2)**2 - c(3)**2 + 2.0d0*c(1)*c(2)*c(3))330:           vol = vol * dsqrt(1.0d0 - c(1)**2 - c(2)**2 - c(3)**2 + 2.0d0*c(1)*c(2)*c(3))
315:        endif331:        endif
316: 332: 
317:        end subroutine get_volume333:        end subroutine get_volume
318: 334: 
319: ! -----------------------------------------------------------------------------------335: ! -----------------------------------------------------------------------------------
320: 336: 
340:       reciplatvec(:,:) = 0.0d0356:       reciplatvec(:,:) = 0.0d0
341:       reciplatvec_grad(:,:,:) = 0.0d0357:       reciplatvec_grad(:,:,:) = 0.0d0
342:       box_lengths(:) = box_params(1:3)358:       box_lengths(:) = box_params(1:3)
343:       box_angles(:) = box_params(4:6)359:       box_angles(:) = box_params(4:6)
344:        360:        
345:       ! cosine of the angles361:       ! cosine of the angles
346:       c(:) = dcos(box_angles(:))362:       c(:) = dcos(box_angles(:))
347:       ! sine of the angles363:       ! sine of the angles
348:       s(:) = dsin(box_angles(:))364:       s(:) = dsin(box_angles(:))
349:       ! factor that is related to the volume (but not quite volume)365:       ! factor that is related to the volume (but not quite volume)
 366:       !print *, 'box_params: ', box_params(:6)
 367:       !print *, 'v: ', v
 368:       !print *, 'v arg: ',(1.0d0 - c(1)**2 - c(2)**2 - c(3)**2 + 2.0d0*c(1)*c(2)*c(3)) 
350:       v = dsqrt(1.0d0 - c(1)**2 - c(2)**2 - c(3)**2 + 2.0d0*c(1)*c(2)*c(3))369:       v = dsqrt(1.0d0 - c(1)**2 - c(2)**2 - c(3)**2 + 2.0d0*c(1)*c(2)*c(3))
351: 370: 
352:       ! define the reciprocal lattice vector matrix371:       ! define the reciprocal lattice vector matrix
353:       ! first row of matrix372:       ! first row of matrix
354:       reciplatvec(1,1) = 1.0d0/box_lengths(1)373:       reciplatvec(1,1) = 1.0d0/box_lengths(1)
355:       ! second row374:       ! second row
356:       reciplatvec(2,1) = -c(3)/(box_lengths(1)*s(3))375:       reciplatvec(2,1) = -c(3)/(box_lengths(1)*s(3))
357:       reciplatvec(2,2) = 1.0d0/(box_lengths(2)*s(3))376:       reciplatvec(2,2) = 1.0d0/(box_lengths(2)*s(3))
358:       ! third row377:       ! third row
359:       reciplatvec(3,1) = (c(3)*(c(1) - c(2)*c(3)) - c(2)*s(3)**2)/(box_lengths(1)*v*s(3))378:       reciplatvec(3,1) = (c(3)*(c(1) - c(2)*c(3)) - c(2)*s(3)**2)/(box_lengths(1)*v*s(3))


r33132/finalio.f90 2017-08-07 15:31:03.225692037 +0100 r33131/finalio.f90 2017-08-07 15:31:06.157731043 +0100
393:            WRITE(MYUNIT2,*) NATOMS393:            WRITE(MYUNIT2,*) NATOMS
394:            WRITE(MYUNIT2, '(2F20.10,I6,1X,E15.8E2)') &394:            WRITE(MYUNIT2, '(2F20.10,I6,1X,E15.8E2)') &
395:                 EDUMMY, LOG_PROD, NSYMOPS, ITDET395:                 EDUMMY, LOG_PROD, NSYMOPS, ITDET
396:            !396:            !
397:         ENDIF397:         ENDIF
398:         !398:         !
399:      ELSE ! <ds656399:      ELSE ! <ds656
400:         WRITE(MYUNIT2,10) J1, QMIN(J1), FF(J1), NPCALL_QMIN(J1)400:         WRITE(MYUNIT2,10) J1, QMIN(J1), FF(J1), NPCALL_QMIN(J1)
401: 10      FORMAT('Energy of minimum ',I6,'=',G20.10, &401: 10      FORMAT('Energy of minimum ',I6,'=',G20.10, &
402:              ' first found at step ',I8,' after ',I20,' function calls')402:              ' first found at step ',I8,' after ',I20,' function calls')
403:         ! dj337: write cell parameters403:         ! dj337
404:         IF (BOXDERIVT) THEN404:         IF (BOXDERIVT) THEN
405:            WRITE(MYUNIT2, *) 'Box lengths: ', boxq(j1,1:3)405:            WRITE(MYUNIT2, *) 'Box lengths: ', boxq(j1,1:3)
406:            IF (.NOT.ORTHO) WRITE(MYUNIT2, *) 'Box angles: ', boxq(j1,4:6)406:            IF (.NOT.ORTHO) WRITE(MYUNIT2, *) 'Box angles: ', boxq(j1,4:6)
407:         ENDIF407:         ENDIF
408:      ENDIF408:      ENDIF
409:      !409:      !
410:      IF (MSORIGT.OR.FRAUSIT) THEN410:      IF (MSORIGT.OR.FRAUSIT) THEN
411:         WRITE(MYUNIT2,20) (QMINP(J1,J2),J2=1,3*(NATOMS-NS))411:         WRITE(MYUNIT2,20) (QMINP(J1,J2),J2=1,3*(NATOMS-NS))
412: 20      FORMAT('Si',3F20.10)412: 20      FORMAT('Si',3F20.10)
413:      ELSE IF (MSTRANST) THEN413:      ELSE IF (MSTRANST) THEN


r33132/finalq.f 2017-08-07 15:31:03.445694967 +0100 r33131/finalq.f 2017-08-07 15:31:06.377733970 +0100
 67:  67: 
 68: ! csw34> Adding in clarification of what is being tightly quenched 68: ! csw34> Adding in clarification of what is being tightly quenched
 69:       WRITE(MYUNIT,'(A)') 'Tightly converging the SAVE lowest energy minima found' 69:       WRITE(MYUNIT,'(A)') 'Tightly converging the SAVE lowest energy minima found'
 70:       WRITE(MYUNIT,'(A)') 'NOTE: these may NOT match the other output files - see below for a sorted list of Lowest minima' 70:       WRITE(MYUNIT,'(A)') 'NOTE: these may NOT match the other output files - see below for a sorted list of Lowest minima'
 71:       DO J1=1,NSAVE 71:       DO J1=1,NSAVE
 72:          IF (QMIN(J1).LT.1.0D10) THEN 72:          IF (QMIN(J1).LT.1.0D10) THEN
 73:             NATOMS=QMINNATOMS(J1) 73:             NATOMS=QMINNATOMS(J1)
 74:             DO J2=1,3*NATOMS 74:             DO J2=1,3*NATOMS
 75:                COORDS(J2,NP)=QMINP(J1,J2) 75:                COORDS(J2,NP)=QMINP(J1,J2)
 76:             ENDDO 76:             ENDDO
 77:             ! dj337: update box parameters with saved ones 
 78:             if (boxderivt) box_params(1:6) = boxq(j1, 1:6) 77:             if (boxderivt) box_params(1:6) = boxq(j1, 1:6)
 79:        78:       
 80:              79:             
 81:             !ds656> 80:             !ds656>
 82:             !write(MYUNIT,*) (QMINT(J1,J2),J2=1,NATOMS) 81:             !write(MYUNIT,*) (QMINT(J1,J2),J2=1,NATOMS)
 83:             IF(NSPECIES(0)>1) THEN 82:             IF(NSPECIES(0)>1) THEN
 84:                CALL SET_ATOMLISTS(QMINT(J1,1:NATOMS),1) 83:                CALL SET_ATOMLISTS(QMINT(J1,1:NATOMS),1)
 85:                CALL SET_LABELS(QMINT(J1,1:NATOMS),NP) 84:                CALL SET_LABELS(QMINT(J1,1:NATOMS),NP)
 86:             ENDIF 85:             ENDIF
 87:             !<ds656 86:             !<ds656
135: !               QMIN(J1)=FULLENERGY134: !               QMIN(J1)=FULLENERGY
136: !            ELSE135: !            ELSE
137: !cl457136: !cl457
138:                QMIN(J1)=POTEL137:                QMIN(J1)=POTEL
139: !            ENDIF138: !            ENDIF
140: 139: 
141:             DO J2=1,3*NATOMS140:             DO J2=1,3*NATOMS
142:                QMINP(J1,J2)=COORDS(J2,NP)141:                QMINP(J1,J2)=COORDS(J2,NP)
143:             ENDDO142:             ENDDO
144: 143: 
145:             ! dj337: save box parameters from end of quench 
146:             if (boxderivt) boxq(j1, 1:6) = box_params(1:6)144:             if (boxderivt) boxq(j1, 1:6) = box_params(1:6)
147: 145: 
148:             IF (CSMT) THEN146:             IF (CSMT) THEN
149:                CSMAV(1:3*NATOMS)=0.0D0147:                CSMAV(1:3*NATOMS)=0.0D0
150:                DO J2=1,CSMGPINDEX148:                DO J2=1,CSMGPINDEX
151: !149: !
152: ! rotate permuted image to best orientation with CSMPMAT150: ! rotate permuted image to best orientation with CSMPMAT
153: ! apply point group operation J2151: ! apply point group operation J2
154: ! 152: ! 
155:                   DO J3=1,NATOMS153:                   DO J3=1,NATOMS


r33132/genrigid.f90 2017-08-07 15:31:03.741698904 +0100 r33131/genrigid.f90 2017-08-07 15:31:06.601736951 +0100
430:   if (boxderivt.and.rigidinit) then430:   if (boxderivt.and.rigidinit) then
431:      xrfrac(1:degfreedoms) = xrigidcoords(1:degfreedoms)431:      xrfrac(1:degfreedoms) = xrigidcoords(1:degfreedoms)
432:      if (ortho) then432:      if (ortho) then
433:         call frac2cart_rb_ortho(nrigidbody, xrigidcoords, xrfrac)433:         call frac2cart_rb_ortho(nrigidbody, xrigidcoords, xrfrac)
434:      else434:      else
435:         call build_H(H, H_grad, .false.)435:         call build_H(H, H_grad, .false.)
436:         call frac2cart_rb_tri(nrigidbody, xrigidcoords, xrfrac, H)436:         call frac2cart_rb_tri(nrigidbody, xrigidcoords, xrfrac, H)
437:      endif437:      endif
438:   endif438:   endif
439: 439: 
 440:   !print *, 'genrigid 437 coords: ', xrigidcoords(1:degfreedoms)
 441: 
440: ! vr274 > are there additional lattice coordinates? If yes, setup transformation matrix442: ! vr274 > are there additional lattice coordinates? If yes, setup transformation matrix
441:   IF(HAS_LATTICE_COORDS) THEN443:   IF(HAS_LATTICE_COORDS) THEN
442:     CALL GET_LATTICE_MATRIX(XRIGIDCOORDS(DEGFREEDOMS-5:DEGFREEDOMS), MLATTICE)444:     CALL GET_LATTICE_MATRIX(XRIGIDCOORDS(DEGFREEDOMS-5:DEGFREEDOMS), MLATTICE)
443:   ELSE ! vr274 > otherwise identity matrix445:   ELSE ! vr274 > otherwise identity matrix
444:     MLATTICE = 0D0446:     MLATTICE = 0D0
445:     MLATTICE(1,1)=1d0447:     MLATTICE(1,1)=1d0
446:     MLATTICE(2,2)=1D0448:     MLATTICE(2,2)=1D0
447:     MLATTICE(3,3)=1D0449:     MLATTICE(3,3)=1D0
448:   ENDIF450:   ENDIF
449: 451: 
595:         ! vr274 > added lattice stuff597:         ! vr274 > added lattice stuff
596:         XRIGIDCOORDS(6*NRIGIDBODY + 3*J1-2:6*NRIGIDBODY + 3*J1) = MATMUL(MLATTICEINV, XCOORDS(3*J9-2:3*J9))598:         XRIGIDCOORDS(6*NRIGIDBODY + 3*J1-2:6*NRIGIDBODY + 3*J1) = MATMUL(MLATTICEINV, XCOORDS(3*J9-2:3*J9))
597:      ENDDO599:      ENDDO
598:   ENDIF600:   ENDIF
599: 601: 
600: ! vr274 > copy lattice coords602: ! vr274 > copy lattice coords
601:   IF(HAS_LATTICE_COORDS) THEN603:   IF(HAS_LATTICE_COORDS) THEN
602:     XRIGIDCOORDS(DEGFREEDOMS - 5:DEGFREEDOMS) =  XCOORDS(3*NATOMS-5:3*NATOMS)604:     XRIGIDCOORDS(DEGFREEDOMS - 5:DEGFREEDOMS) =  XCOORDS(3*NATOMS-5:3*NATOMS)
603:   ENDIF605:   ENDIF
604: 606: 
 607:   !print *, 'transformctorigid rb coords'
 608:   !do j1 = 1, 2*nrigidbody
 609:   !   print *, xrigidcoords(3*j1-2:3*j1)
 610:   !enddo
 611: 
605:   ! dj337: if computing box derivatives, convert rb coords to fractional612:   ! dj337: if computing box derivatives, convert rb coords to fractional
606:   if (boxderivt.and.rigidinit) then613:   if (boxderivt.and.rigidinit) then
607:      if (ortho) then614:      if (ortho) then
608:         call cart2frac_rb_ortho(nrigidbody, xrigidcoords, xrfrac)615:         call cart2frac_rb_ortho(nrigidbody, xrigidcoords, xrfrac)
609:      else616:      else
610:         call build_H(H, H_grad, .false.)617:         call build_H(H, H_grad, .false.)
611:         call inversematrix(H, H_inv)618:         call inversematrix(H, H_inv)
612:         call cart2frac_rb_tri(nrigidbody, xrigidcoords, xrfrac, H_inv)619:         call cart2frac_rb_tri(nrigidbody, xrigidcoords, xrfrac, H_inv)
613:      endif620:      endif
614:      xrigidcoords(1:degfreedoms) = xrfrac(1:degfreedoms)621:      xrigidcoords(1:degfreedoms) = xrfrac(1:degfreedoms)
615:   endif622:   endif
616: 623: 
 624:   !print *, 'transformctorigid rb frac coords'
 625:   !do j1 = 1, 2*nrigidbody
 626:   !   print *, xrigidcoords(3*j1-2:3*j1)
 627:   !enddo 
 628: 
617: END SUBROUTINE TRANSFORMCTORIGID629: END SUBROUTINE TRANSFORMCTORIGID
618: 630: 
619: !-----------------------------------------------------------631: !-----------------------------------------------------------
620: 632: 
621: SUBROUTINE TRANSFORMCTORIGID_OLD (XCOORDS, XRIGIDCOORDS)633: SUBROUTINE TRANSFORMCTORIGID_OLD (XCOORDS, XRIGIDCOORDS)
622: 634: 
623:   USE COMMONS, ONLY: NATOMS635:   USE COMMONS, ONLY: NATOMS
624:   USE VEC3636:   USE VEC3
625:   IMPLICIT NONE637:   IMPLICIT NONE
626: 638: 


r33132/keywords.f 2017-08-07 15:31:04.053703054 +0100 r33131/keywords.f 2017-08-07 15:31:06.833740034 +0100
3623:          ENDIF3623:          ENDIF
3624: 3624: 
3625:          BENZRIGIDEWALDT = .TRUE.3625:          BENZRIGIDEWALDT = .TRUE.
3626: 3626: 
3627:          IF (NITEMS.GT.1) CALL READF(EWALDREALC)3627:          IF (NITEMS.GT.1) CALL READF(EWALDREALC)
3628:          IF (NITEMS.GT.2) CALL READF(EWALDRECIPC)3628:          IF (NITEMS.GT.2) CALL READF(EWALDRECIPC)
3629: 3629: 
3630:          ! CALCULATE ALPHA = 5.6/L_MIN3630:          ! CALCULATE ALPHA = 5.6/L_MIN
3631:          EWALDALPHA = 5.6D0/12.0d03631:          EWALDALPHA = 5.6D0/12.0d0
3632: 3632: 
 3633:          ! SET NUMBER OF LATTICE AND RECIPROCAL LATTICE VECTORS
 3634:          !NEWALDREAL(1) = FLOOR(EWALDREALC/BOXLX+0.5D0)
 3635:          !NEWALDREAL(2) = FLOOR(EWALDREALC/BOXLY+0.5D0)
 3636:          !NEWALDREAL(3) = FLOOR(EWALDREALC/BOXLZ+0.5D0)
 3637: 
 3638:          !NEWALDRECIP(1) = FLOOR(EWALDRECIPC*BOXLX/(2*PI))
 3639:          !NEWALDRECIP(2) = FLOOR(EWALDRECIPC*BOXLY/(2*PI))
 3640:          !NEWALDRECIP(3) = FLOOR(EWALDRECIPC*BOXLZ/(2*PI))
 3641: 
 3642:          ! ALLOCATE ARRAYS FOR STRUCTURE FACTORS
 3643:          !ALLOCATE(RERHOARRAY(2*NEWALDRECIP(1)+1, 2*NEWALDRECIP(2)+1, 2*NEWALDRECIP(3)+1))
 3644:          !ALLOCATE(IMRHOARRAY(2*NEWALDRECIP(1)+1, 2*NEWALDRECIP(2)+1, 2*NEWALDRECIP(3)+1))
 3645: 
3633:          ! ALLOCATE BENZENE MOLECULE3646:          ! ALLOCATE BENZENE MOLECULE
3634:          NRBSITES = 123647:          NRBSITES = 12
3635:          ALLOCATE(SITE(NRBSITES,3))3648:          ALLOCATE(SITE(NRBSITES,3))
3636:          ALLOCATE(RBSTLA(NRBSITES,3))3649:          ALLOCATE(RBSTLA(NRBSITES,3))
3637:          CALL DEFPAHARIGID()3650:          CALL DEFPAHARIGID()
3638:          NCARBON = 63651:          NCARBON = 6
3639: 3652: 
3640: ! ----------------------------------------------------------------------------------------3653: ! ----------------------------------------------------------------------------------------
3641:       ! dj337: Ewald summation3654:       ! dj337: Ewald summation
3642:       ELSE IF (WORD.EQ.'EWALD') THEN3655:       ELSE IF (WORD.EQ.'EWALD') THEN
3646:             STOP3659:             STOP
3647:          ENDIF3660:          ENDIF
3648: 3661: 
3649:          EWALDT = .TRUE.3662:          EWALDT = .TRUE.
3650: 3663: 
3651:          IF (NITEMS.GT.1) CALL READI(EWALDN)3664:          IF (NITEMS.GT.1) CALL READI(EWALDN)
3652:          IF (NITEMS.GT.2) CALL READF(EWALDREALC)3665:          IF (NITEMS.GT.2) CALL READF(EWALDREALC)
3653:          IF (NITEMS.GT.3) CALL READF(EWALDRECIPC)3666:          IF (NITEMS.GT.3) CALL READF(EWALDRECIPC)
3654:          IF (NITEMS.GT.4) CALL READF(RSPEED)3667:          IF (NITEMS.GT.4) CALL READF(RSPEED)
3655: 3668: 
 3669:          !VOL = BOXLX*BOXLY*BOXLZ
 3670:          !CALL VOLUME(VOL)
 3671:          ! calculate alpha value
 3672: !        ! EWALDALPHA = 0.28D0
 3673:          !EWALDALPHA = (RSPEED*(PI**3)*NATOMS/(VOL**2))**(1.0D0/6.0D0)
 3674: 
3656:          ewaldalpha = 5.6d0/minval(box3d)3675:          ewaldalpha = 5.6d0/minval(box3d)
3657: 3676: 
 3677:          !IF (ORTHO) THEN
 3678:          !   ! set number of lattice vectors
 3679:          !   NEWALDREAL(1) = FLOOR(EWALDREALC/BOXLX+0.5D0)
 3680:          !   NEWALDREAL(2) = FLOOR(EWALDREALC/BOXLY+0.5D0)
 3681:          !   NEWALDREAL(3) = FLOOR(EWALDREALC/BOXLZ+0.5D0)
 3682: 
 3683:          !   ! set number of reciprocal lattice vectors
 3684:          !   NEWALDRECIP(1) = FLOOR(EWALDRECIPC*BOXLX/(2*PI))
 3685:          !   NEWALDRECIP(2) = FLOOR(EWALDRECIPC*BOXLY/(2*PI))
 3686:          !   NEWALDRECIP(3) = FLOOR(EWALDRECIPC*BOXLZ/(2*PI))
 3687:          !ELSE
 3688:          !   print *, 'Not yet implemented for non-orthorhombic boxes!'
 3689:          !   STOP
 3690:          !ENDIF
 3691: 
 3692:          !ALLOCATE(RERHOARRAY(2*NEWALDRECIP(1)+1, 2*NEWALDRECIP(2)+1, 2*NEWALDRECIP(3)+1))
 3693:          !ALLOCATE(IMRHOARRAY(2*NEWALDRECIP(1)+1, 2*NEWALDRECIP(2)+1, 2*NEWALDRECIP(3)+1))
 3694: 
 3695: 
3658:          NRBSITES = 123696:          NRBSITES = 12
3659:          ALLOCATE(SITE(NRBSITES,3))3697:          ALLOCATE(SITE(NRBSITES,3))
3660:          ALLOCATE(RBSTLA(NRBSITES,3))3698:          ALLOCATE(RBSTLA(NRBSITES,3))
 3699:          !ALLOCATE(STCHRG(2*NRBSITES))
3661: 3700: 
3662:          CALL DEFPAHARIGID()3701:          CALL DEFPAHARIGID()
3663:          NCARBON = 63702:          NCARBON = 6
3664:          !nrigidbody = 323703:          nrigidbody = 32
3665:          !CALL DEFBENZENERIGID()3704:          !CALL DEFBENZENERIGID()
3666: 3705: 
3667: ! ----------------------------------------------------------------------------------------3706: ! ----------------------------------------------------------------------------------------
3668: 3707: 
3669:       ELSE IF (WORD.EQ.'FAL') THEN3708:       ELSE IF (WORD.EQ.'FAL') THEN
3670:          FAL=.TRUE.3709:          FAL=.TRUE.
3671: 3710: 
3672:       ELSE IF (WORD == 'FEBH') THEN3711:       ELSE IF (WORD == 'FEBH') THEN
3673:          CALL READF(FETEMP)3712:          CALL READF(FETEMP)
3674:          FEBHT = .TRUE.3713:          FEBHT = .TRUE.


r33132/main.F 2017-08-07 15:31:04.313706512 +0100 r33131/main.F 2017-08-07 15:31:07.053742964 +0100
194:           CALL GENRIGID_READ_FROM_FILE ()194:           CALL GENRIGID_READ_FROM_FILE ()
195: ! csw34> Tell the user how many degrees of freedom there are in the system195: ! csw34> Tell the user how many degrees of freedom there are in the system
196:           WRITE(MYUNIT, '(A,I15)') " genrigid> rbodyconfig used to specifiy rigid bodies, degrees of freedom now ", DEGFREEDOMS196:           WRITE(MYUNIT, '(A,I15)') " genrigid> rbodyconfig used to specifiy rigid bodies, degrees of freedom now ", DEGFREEDOMS
197:           IF (GCBHT) THEN197:           IF (GCBHT) THEN
198: ! csw34> Make sure we aren't running GCBH with rigid bodies - very dangerous! Could be done but only with great care...198: ! csw34> Make sure we aren't running GCBH with rigid bodies - very dangerous! Could be done but only with great care...
199:               WRITE(MYUNIT, '(A)') " genrigid> ERROR: cannot use rigid bodies with GCBH! Stopping."199:               WRITE(MYUNIT, '(A)') " genrigid> ERROR: cannot use rigid bodies with GCBH! Stopping."
200:               STOP                200:               STOP                
201:           END IF201:           END IF
202:       END IF202:       END IF
203: 203: 
 204:       !print *, sitesrigidbody(1, 1:3, 1)
 205:       !print *, 'coords: ', size(coords,1), size(coords,2)
 206: 
204: ! dj337: allocate charges for benzenes207: ! dj337: allocate charges for benzenes
205:       if (benzrigidewaldt.or.ewaldt) then208:       if (benzrigidewaldt.or.ewaldt) then
206:          ALLOCATE(STCHRG(NRIGIDBODY*NRBSITES))209:          ALLOCATE(STCHRG(NRIGIDBODY*NRBSITES))
207:          CALL DEFBENZENERIGIDEWALD()210:          CALL DEFBENZENERIGIDEWALD()
208:       endif211:       endif
209: 212: 
210: 213: 
211:       IF (MULTIPOTT) THEN214:       IF (MULTIPOTT) THEN
212:           CALL MULTIPOT_INITIALISE()215:           CALL MULTIPOT_INITIALISE()
213:       ENDIF216:       ENDIF
344: !        though - but only from xyz format! Is there a hidden347: !        though - but only from xyz format! Is there a hidden
345: !        conversion?348: !        conversion?
346: !         349: !         
347:          OPEN(UNIT=1,FILE='compare',STATUS='OLD')350:          OPEN(UNIT=1,FILE='compare',STATUS='OLD')
348:          READ(1,*) (COORCOMP(J1),J1=1,3*NATOMS)351:          READ(1,*) (COORCOMP(J1),J1=1,3*NATOMS)
349:          CLOSE(1)352:          CLOSE(1)
350:       ENDIF353:       ENDIF
351: 354: 
352:       ALLOCATE(FF(NSAVE),QMIN(MAX(NSAVE,1)))355:       ALLOCATE(FF(NSAVE),QMIN(MAX(NSAVE,1)))
353:       ALLOCATE(QMINP(NSAVE,3*NATOMSALLOC))356:       ALLOCATE(QMINP(NSAVE,3*NATOMSALLOC))
354:       ! dj337: allocate array for saving box parameters 
355:       if (boxderivt) allocate(boxq(nsave,6))357:       if (boxderivt) allocate(boxq(nsave,6))
356:       ALLOCATE(QMINNATOMS(NSAVE))358:       ALLOCATE(QMINNATOMS(NSAVE))
357:       ALLOCATE(AVOIDNATOMS(MAXSAVE))359:       ALLOCATE(AVOIDNATOMS(MAXSAVE))
358:       ALLOCATE(QMINT(NSAVE,NATOMSALLOC))360:       ALLOCATE(QMINT(NSAVE,NATOMSALLOC))
359:       ALLOCATE(NPCALL_QMIN(NSAVE))361:       ALLOCATE(NPCALL_QMIN(NSAVE))
360:       IF (MONITORT) THEN362:       IF (MONITORT) THEN
361:          ALLOCATE(LOWESTC(3*NATOMSALLOC))363:          ALLOCATE(LOWESTC(3*NATOMSALLOC))
362:       ENDIF364:       ENDIF
363: 365: 
364: !        csw34> ALLOCATE the interaction energy tracking arrays if A9INTE in data366: !        csw34> ALLOCATE the interaction energy tracking arrays if A9INTE in data


r33132/mc.F 2017-08-07 15:31:04.573709968 +0100 r33131/mc.F 2017-08-07 15:31:07.277745940 +0100
482: ! jwrm2> end482: ! jwrm2> end
483:          EPPREV(JP)=0.0D0483:          EPPREV(JP)=0.0D0
484:          ELASTSYM(JP)=0.0D0484:          ELASTSYM(JP)=0.0D0
485:          IF (.NOT.RESTORET) EBEST(JP)=POTEL485:          IF (.NOT.RESTORET) EBEST(JP)=POTEL
486:          BESTCOORDS(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)486:          BESTCOORDS(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)
487:          JBEST(JP)=0487:          JBEST(JP)=0
488:          RMIN=POTEL488:          RMIN=POTEL
489:          RCOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,1)489:          RCOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,1)
490:          COORDSO(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)490:          COORDSO(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)
491:          if (boxderivt) box_paramso(1:6) = box_params(1:6) ! dj337: save box parameters491:          if (boxderivt) box_paramso(1:6) = box_params(1:6) ! dj337: save box parameters
 492:          !print *, 'intl saving'
 493:          !print *, 'boxparamso: ', box_paramso(1:6)
 494:          !if (boxderivt) tempbox_params(1:6) = box_params(1:6)
 495:          ! print *, 'coordso   : ', coordso(1:25, jp)
492:          LABELSO(1:NATOMS,JP)=LABELS(1:NATOMS,JP) ! <ds656496:          LABELSO(1:NATOMS,JP)=LABELS(1:NATOMS,JP) ! <ds656
493:          VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)497:          VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)
494: !hk286 - otherwise NEWRESTART is always called when the keyword RESTORE is also used498: !hk286 - otherwise NEWRESTART is always called when the keyword RESTORE is also used
495:          IF (.NOT.RESTORET) THEN499:          IF (.NOT.RESTORET) THEN
496:             EBEST(JP)=POTEL500:             EBEST(JP)=POTEL
497:             JBEST(JP)=0501:             JBEST(JP)=0
498:             BESTCOORDS(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)502:             BESTCOORDS(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)
499:          ENDIF503:          ENDIF
500: !504: !
501: ! Set up data for relxation blocks with grand canonical BH505: ! Set up data for relxation blocks with grand canonical BH
827:                ENDIF831:                ENDIF
828: 832: 
829:                IF(ROTATEHINGET .AND. MOD(J1,ROTATEHINGEFREQ).EQ.0) THEN833:                IF(ROTATEHINGET .AND. MOD(J1,ROTATEHINGEFREQ).EQ.0) THEN
830:                        DOROTATEHINGE=.TRUE.834:                        DOROTATEHINGE=.TRUE.
831:                ENDIF835:                ENDIF
832: 836: 
833: !837: !
834: ! csw34> Coordinates are saved so that moves can be undone838: ! csw34> Coordinates are saved so that moves can be undone
835: !839: !
836:                SAVECOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,JP)840:                SAVECOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,JP)
837:                if (boxderivt) savebox_params(1:6) = box_params(1:6) ! dj337: save box params841:                if (boxderivt) savebox_params(1:6) = box_params(1:6) ! save box params
838: ! csw34> If you want to look at the effect of moves, you can dump out842: ! csw34> If you want to look at the effect of moves, you can dump out
839: ! the structure BEFORE the move here.843: ! the structure BEFORE the move here.
840: !                 CALL A9DUMPPDB(COORDS(:,JP),"beforemove")844: !                 CALL A9DUMPPDB(COORDS(:,JP),"beforemove")
841: !                 CALL CHARMMDUMP(COORDS(:,JP),'beforemove')845: !                 CALL CHARMMDUMP(COORDS(:,JP),'beforemove')
842:                IF (MACROIONT) THEN846:                IF (MACROIONT) THEN
843:                   IF(RANDOMSEEDT) THEN847:                   IF(RANDOMSEEDT) THEN
844:                     CALL DATE_AND_TIME(DATECHAR,TIMECHAR,ZONECHAR,VALUES)848:                     CALL DATE_AND_TIME(DATECHAR,TIMECHAR,ZONECHAR,VALUES)
845:                     ITIME1 = VALUES(6) * 60 + VALUES(7)849:                     ITIME1 = VALUES(6) * 60 + VALUES(7)
846:                     CALL SDPRND(ITIME1 + MYNODE)850:                     CALL SDPRND(ITIME1 + MYNODE)
847:                   END IF851:                   END IF
1462:                   ENDIF1466:                   ENDIF
1463: ! END OF KEYWORD <DUMPSTEPS> BLOCK1467: ! END OF KEYWORD <DUMPSTEPS> BLOCK
1464: 1468: 
1465: 1469: 
1466:                NQ(JP)=NQ(JP)+11470:                NQ(JP)=NQ(JP)+1
1467:                IF(CHRMMT.AND.ACESOLV) NCHENCALLS=ACEUPSTEP-11471:                IF(CHRMMT.AND.ACESOLV) NCHENCALLS=ACEUPSTEP-1
1468: !fh301>{{{1472: !fh301>{{{
1469:                STRUC = NQ(JP)1473:                STRUC = NQ(JP)
1470: !fh301>}}}1474: !fh301>}}}
1471:                CALL QUENCH(.FALSE.,JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)  1475:                CALL QUENCH(.FALSE.,JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)  
1472:                ! dj337: save cell parameters from end of quench in a temp array 
1473:                if (boxderivt) tempbox_params(1:6) = box_params(1:6)1476:                if (boxderivt) tempbox_params(1:6) = box_params(1:6)
 1477:                !print *, 'potel, eprev: ', potel, eprev(jp)
 1478:                !print *, 'boxparams: ', box_params(1:6)
 1479:                !if (boxderivt.and.atest) then
 1480:                !   print *, 'potel, eprev: ', potel, eprev(jp)
 1481:                !   tempbox_params(1:6) = box_params(1:6) ! dj337: save box params
 1482:                !   print *, 'tempbox: ', tempbox_params(1:6)
 1483:                !endif
 1484:                !print *, 'energy post-quench: ', eprev(jp), potel
 1485:                !print *, 'saving params'
 1486:                !print *, 'tempbox: ', tempbox_params(1:6)
1474:                NQTOT=NQTOT+11487:                NQTOT=NQTOT+1
1475: 1488: 
1476: !  Check for results of taboo list. SELFT is set in taboo.1489: !  Check for results of taboo list. SELFT is set in taboo.
1477:  1490:  
1478:                IF (SELFT) THEN1491:                IF (SELFT) THEN
1479:                   CALL MYRESET(JP,NATOMS,NPAR,NSEED)1492:                   CALL MYRESET(JP,NATOMS,NPAR,NSEED)
1480:                   IF (DEBUG) THEN1493:                   IF (DEBUG) THEN
1481:                      WRITE(MYUNIT,'(A)') 'Taboo list:'1494:                      WRITE(MYUNIT,'(A)') 'Taboo list:'
1482:                      WRITE(MYUNIT,'(6F20.10)') (ESAVE(J2,1),J2=1,NT(1))1495:                      WRITE(MYUNIT,'(6F20.10)') (ESAVE(J2,1),J2=1,NT(1))
1483:                      WRITE(MYUNIT,73) JP,J1,POTEL,EPREV(JP),NSUCCESS(JP),NFAIL(JP)1496:                      WRITE(MYUNIT,73) JP,J1,POTEL,EPREV(JP),NSUCCESS(JP),NFAIL(JP)
1519:      1                       ' Markov E=',EPREV(JP),' t=',TIME-TSTART1532:      1                       ' Markov E=',EPREV(JP),' t=',TIME-TSTART
1520:                      ENDIF1533:                      ENDIF
1521:                   ELSE1534:                   ELSE
1522: !     fh301>}}}1535: !     fh301>}}}
1523:                      IF (NPAR.GT.1) THEN1536:                      IF (NPAR.GT.1) THEN
1524:                         IF (.NOT.SUPPRESST) WRITE(MYUNIT,'(A,I0.2,A,I10,A,G20.10,A,I5,A,G12.5,A,G20.10,A,F11.1)') '[',JP,1537:                         IF (.NOT.SUPPRESST) WRITE(MYUNIT,'(A,I0.2,A,I10,A,G20.10,A,I5,A,G12.5,A,G20.10,A,F11.1)') '[',JP,
1525:      1                       ']Qu ',NQ(JP),' E=',POTEL,' steps=',ITERATIONS,' RMS=',RMS,' Markov E=',EPREV(JP),' t=',TIME-TSTART1538:      1                       ']Qu ',NQ(JP),' E=',POTEL,' steps=',ITERATIONS,' RMS=',RMS,' Markov E=',EPREV(JP),' t=',TIME-TSTART
1526:                      ELSE1539:                      ELSE
1527:                         IF (.NOT.SUPPRESST) WRITE(MYUNIT,'(A,I10,A,G20.10,A,I5,A,G12.5,A,G20.10,A,F11.1)') 'Qu ',NQ(JP),' E=',1540:                         IF (.NOT.SUPPRESST) WRITE(MYUNIT,'(A,I10,A,G20.10,A,I5,A,G12.5,A,G20.10,A,F11.1)') 'Qu ',NQ(JP),' E=',
1528:      1                       POTEL,' steps=',ITERATIONS,' RMS=',RMS,' Markov E=',EPREV(JP),' t=',TIME-TSTART1541:      1                       POTEL,' steps=',ITERATIONS,' RMS=',RMS,' Markov E=',EPREV(JP),' t=',TIME-TSTART
1529:                         ! dj337: write box parameters to output file 
1530:                         if (boxderivt) write(myunit, *) 'Box parameters:1542:                         if (boxderivt) write(myunit, *) 'Box parameters:
1531:      1                     ', box_params(1:6)1543:      1                     ', box_params(1:6)
1532: !fh301>{{{1544: !fh301>{{{
1533:                      ENDIF1545:                      ENDIF
1534: !fh301>}}}1546: !fh301>}}}
1535:                   ENDIF1547:                   ENDIF
1536:                   1548:                   
1537: !                       1549: !                       
1538: ! khs26> Printing energy decomposition for AMBER energies if we're using AMBER.1550: ! khs26> Printing energy decomposition for AMBER energies if we're using AMBER.
1539:                   IF (ENERGY_DECOMPT) THEN1551:                   IF (ENERGY_DECOMPT) THEN
1991: !2003: !
1992: !  Sanity check to make sure the Markov energy agrees with COORDSO. 2004: !  Sanity check to make sure the Markov energy agrees with COORDSO. 
1993: !  Stop if not true.2005: !  Stop if not true.
1994: !2006: !
1995:                IF ((DEBUG.OR.CHECKMARKOVT).AND.(.NOT.(HBONDMATRIX.OR.PERMOPT.OR.PERMINVOPT.OR.FEBHT.OR.QALCST.OR. 2007:                IF ((DEBUG.OR.CHECKMARKOVT).AND.(.NOT.(HBONDMATRIX.OR.PERMOPT.OR.PERMINVOPT.OR.FEBHT.OR.QALCST.OR. 
1996:      &                           MODULART.OR.PERCGROUPT))) THEN2008:      &                           MODULART.OR.PERCGROUPT))) THEN
1997: ! jwrm2> If we're using percolate, we need to set compression to whether it was on or not when calculating the Markov energy2009: ! jwrm2> If we're using percolate, we need to set compression to whether it was on or not when calculating the Markov energy
1998:                   QUENCHCOMP = COMPON2010:                   QUENCHCOMP = COMPON
1999:                   IF (PERCOLATET) COMPON = PERCCOMPMARKOV2011:                   IF (PERCOLATET) COMPON = PERCCOMPMARKOV
2000:                   if (boxderivt) box_params(1:6) = box_paramso(1:6) ! dj337: restore box parameters2012:                   if (boxderivt) box_params(1:6) = box_paramso(1:6) ! dj337: restore box parameters
 2013:                   !print *, 'calling OPOTEL!'
 2014:                   !print *, 'boxparamso: ', box_paramso(1:6)
 2015:                   !print *, 'coordso   : ', coordso(1:25, jp)
2001:                   CALL POTENTIAL(COORDSO(:,JP),GRAD,OPOTEL,.FALSE.,.FALSE.)2016:                   CALL POTENTIAL(COORDSO(:,JP),GRAD,OPOTEL,.FALSE.,.FALSE.)
2002:                   if (boxderivt) box_params(1:6) = tempbox_params(1:6) ! dj337: put them back2017:                   if (boxderivt) box_params(1:6) = tempbox_params(1:6)
2003:                   IF (ABS(OPOTEL-EPREV(JP)).GT.ABS(ECONV)) THEN2018:                   IF (ABS(OPOTEL-EPREV(JP)).GT.ABS(ECONV)) THEN
2004:                      IF (EVAP) THEN2019:                      IF (EVAP) THEN
2005:                         WRITE(MYUNIT,'(3(A,G20.10))') 'mc> WARNING - energy for saved coordinates ',OPOTEL,2020:                         WRITE(MYUNIT,'(3(A,G20.10))') 'mc> WARNING - energy for saved coordinates ',OPOTEL,
2006:      &                     ' differs from Markov energy ',EPREV(JP),' because an atom moved outside the container'2021:      &                     ' differs from Markov energy ',EPREV(JP),' because an atom moved outside the container'
2007:                      ELSE2022:                      ELSE
2008:                         WRITE(MYUNIT,'(2(A,G20.10))') 'mc> ERROR - energy for coordinates in COORDSO=',OPOTEL,2023:                         WRITE(MYUNIT,'(2(A,G20.10))') 'mc> ERROR - energy for coordinates in COORDSO=',OPOTEL,
2009:      &                                                 ' but Markov energy=',EPREV(JP)2024:      &                                                 ' but Markov energy=',EPREV(JP)
2010:                         IF (.NOT.(CHEMSHIFT2)) STOP2025:                         IF (.NOT.(CHEMSHIFT2)) STOP
2011:                      ENDIF2026:                      ENDIF
2012:                   ENDIF2027:                   ENDIF
2013:                ENDIF2028:                ENDIF
2014: 2029: 
 2030:                ! dj337
 2031:                !print *, 'giving saved params'
 2032:                !if (boxderivt) box_paramso(1:6) = tempbox_params(1:6)
 2033:                !print *, 'boxparamso: ', box_paramso(1:6)
 2034:                
2015:                IF (RATIOT) THEN2035:                IF (RATIOT) THEN
2016:                   CALL NULLMOVE(SCREENC,COORDSO(:,JP),ATEST,NULLMOVES,JP)2036:                   CALL NULLMOVE(SCREENC,COORDSO(:,JP),ATEST,NULLMOVES,JP)
2017:                ENDIF2037:                ENDIF
2018: ! Accept or reject step. If the quench did not converge then allow a2038: ! Accept or reject step. If the quench did not converge then allow a
2019: ! potenial move, but count it as a rejection in terms of NSUCCESS and2039: ! potenial move, but count it as a rejection in terms of NSUCCESS and
2020: ! NFAIL. This way we will accept a lower minimum if found, but the steps won;t become so big.2040: ! NFAIL. This way we will accept a lower minimum if found, but the steps won;t become so big.
2021: ! However, for TIP5P some cold fusion events that had not actually reached the threshold for2041: ! However, for TIP5P some cold fusion events that had not actually reached the threshold for
2022: ! rejection were actually accepted. Must prevent this!2042: ! rejection were actually accepted. Must prevent this!
2023:                IF (ATEST) THEN2043:                IF (ATEST) THEN
2024: #ifdef MPI2044: #ifdef MPI
2080:                            HBONDGROUPS(HBONDMATCHN,:,:)=HBONDMAT2100:                            HBONDGROUPS(HBONDMATCHN,:,:)=HBONDMAT
2081: ! csw34> Dump an updates group*.mat file for the group2101: ! csw34> Dump an updates group*.mat file for the group
2082:                            WRITE(NHBONDGROUPSCHAR,*) HBONDMATCHN2102:                            WRITE(NHBONDGROUPSCHAR,*) HBONDMATCHN
2083:                            HBONDGROUPNAME='group'//TRIM(ADJUSTL(NHBONDGROUPSCHAR))2103:                            HBONDGROUPNAME='group'//TRIM(ADJUSTL(NHBONDGROUPSCHAR))
2084:                            CALL WRITEHBONDMATRIX(HBONDMAT(:,:),HBONDNRES,HBONDGROUPNAME)2104:                            CALL WRITEHBONDMATRIX(HBONDMAT(:,:),HBONDNRES,HBONDGROUPNAME)
2085:                         ENDIF2105:                         ENDIF
2086:                      ENDIF2106:                      ENDIF
2087:                   ENDIF2107:                   ENDIF
2088: ! END OF HBONDMATRIX2108: ! END OF HBONDMATRIX
2089:                   COORDSO(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)2109:                   COORDSO(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)
2090:                   if (boxderivt) box_paramso(1:6) = box_params(1:6) ! dj337: update saved box parameters2110:                   if (boxderivt) box_paramso(1:6) = box_params(1:6)
 2111:                   !if (boxderivt) box_paramso(1:6) = box_params(1:6) ! dj337: update saved box params
 2112:                   !print *, 'boxparamso: ', box_paramso(1:6)
 2113:                   !print *, 'coordso   : ', coordso(1:25, jp)
2091:                   LABELSO(1:NATOMS,JP)=LABELS(1:NATOMS,JP) ! <ds6562114:                   LABELSO(1:NATOMS,JP)=LABELS(1:NATOMS,JP) ! <ds656
2092:                   VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)2115:                   VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)
2093:                ELSE2116:                ELSE
2094:                   NFAIL(JP)=NFAIL(JP)+12117:                   NFAIL(JP)=NFAIL(JP)+1
2095:                   REJSTREAK=REJSTREAK+12118:                   REJSTREAK=REJSTREAK+1
2096:                   MARKOVSTEP=.FALSE.2119:                   MARKOVSTEP=.FALSE.
2097:                   CALL MYRESET(JP,NATOMS,NPAR,NSEED)2120:                   CALL MYRESET(JP,NATOMS,NPAR,NSEED)
2098:                   IF (DEBUG) THEN2121:                   IF (DEBUG) THEN
2099:                      WRITE(MYUNIT,36) JP,RANDOM,POTEL,EPREV(JP),NSUCCESS(JP),NFAIL(JP)2122:                      WRITE(MYUNIT,36) JP,RANDOM,POTEL,EPREV(JP),NSUCCESS(JP),NFAIL(JP)
2100: 36                   FORMAT('JP,RAN,POTEL,EPREV,NSUC,NFAIL=',I2,3F15.7,2I6,' REJ')2123: 36                   FORMAT('JP,RAN,POTEL,EPREV,NSUC,NFAIL=',I2,3F15.7,2I6,' REJ')


r33132/potential.f90 2017-08-07 15:31:04.821713271 +0100 r33131/potential.f90 2017-08-07 15:31:07.497748868 +0100
126: !   END IF126: !   END IF
127: 127: 
128: ! ds656> For PTMC runs with MIEF and BOXCENTROID, safer to128: ! ds656> For PTMC runs with MIEF and BOXCENTROID, safer to
129: !        call BOXCENTROID on every potential call...129: !        call BOXCENTROID on every potential call...
130:    IF(BOXCENTROIDT .AND. PTMC) CALL BOXCENTROID(X)130:    IF(BOXCENTROIDT .AND. PTMC) CALL BOXCENTROID(X)
131: !131: !
132: ! Count the number of potential calls.132: ! Count the number of potential calls.
133:    NPCALL=NPCALL+1133:    NPCALL=NPCALL+1
134: 10 CONTINUE134: 10 CONTINUE
135: 135: 
 136:    !print *, 'intl coords pot: '
 137:    !do j1 = 1, 2*nrigidbody
 138:    !   print *, x(3*j1-2:3*j1)
 139:    !enddo
 140: 
136:    ! dj337: box derivatives141:    ! dj337: box derivatives
137:    if (boxderivt) then142:    if (boxderivt) then
138:       ! store updated cell parameters143:       ! store updated cell parameters
139:       boxlx = box_params(1)144:       boxlx = box_params(1)
140:       boxly = box_params(2)145:       boxly = box_params(2)
141:       boxlz = box_params(3)146:       boxlz = box_params(3)
142:       palpha = box_params(4)147:       palpha = box_params(4)
143:       pbeta = box_params(5)148:       pbeta = box_params(5)
144:       pgamma = box_params(6)149:       pgamma = box_params(6)
145:       ! reset lattice derivatives150:       ! reset lattice derivatives
147:       ! convert fractional to absolute coordinates152:       ! convert fractional to absolute coordinates
148:       ! NOTE: don't convert if rigid bodies and already153:       ! NOTE: don't convert if rigid bodies and already
149:       ! in atomic coordinates, cause those are absolute154:       ! in atomic coordinates, cause those are absolute
150:       ! by default when rigidinit155:       ! by default when rigidinit
151:       if (rigidinit.and.(.not.atomrigidcoordt)) then156:       if (rigidinit.and.(.not.atomrigidcoordt)) then
152:          xfrac(1:3*natoms) = x(1:3*natoms)157:          xfrac(1:3*natoms) = x(1:3*natoms)
153:          call frac2cart(x, xfrac)158:          call frac2cart(x, xfrac)
154:       endif159:       endif
155:    endif 160:    endif 
156: 161: 
 162:    !print *, 'final coords pot: '
 163:    !do j1 = 1, 2*nrigidbody
 164:    !   print *, x(3*j1-2:3*j1)
 165:    !enddo
 166: 
157:    IF (MULTIPOTT) THEN167:    IF (MULTIPOTT) THEN
158:        IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN168:        IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
159:            XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)169:            XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)
160:            CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)170:            CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)
161:        END IF171:        END IF
162: !       CALL RAD(X, GRAD, EREAL, GRADT)172: !       CALL RAD(X, GRAD, EREAL, GRADT)
163:        CALL MULTIPOT_CALL(X, GRADATOMS, EREAL, GRADT, SECT)173:        CALL MULTIPOT_CALL(X, GRADATOMS, EREAL, GRADT, SECT)
164:        GRAD(1:3*NATOMS)=GRADATOMS(:)174:        GRAD(1:3*NATOMS)=GRADATOMS(:)
165: 175: 
166:        IF ( RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN176:        IF ( RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
1039:       CALL RAD(X, GRAD, EREAL, GRADT)1049:       CALL RAD(X, GRAD, EREAL, GRADT)
1040:       CALL MLJ(X, GRAD, EREAL, GRADT, SECT, .FALSE.)1050:       CALL MLJ(X, GRAD, EREAL, GRADT, SECT, .FALSE.)
1041: 1051: 
1042:    ELSE IF (GLJT) THEN1052:    ELSE IF (GLJT) THEN
1043: ! generalised LJ with no cutoff1053: ! generalised LJ with no cutoff
1044:       CALL RAD(X, GRAD, EREAL, GRADT)1054:       CALL RAD(X, GRAD, EREAL, GRADT)
1045:       CALL GLJ(X, GRAD, EREAL, GRADT)1055:       CALL GLJ(X, GRAD, EREAL, GRADT)
1046: 1056: 
1047:    ELSE IF (BINARY) THEN1057:    ELSE IF (BINARY) THEN
1048:       IF (SHIFTCUT) THEN1058:       IF (SHIFTCUT) THEN
1049:          ! dj377: shifcut without minimum image convention1059:          !print *, 'in shiftcut!'
1050:          if (rigidinit.and.(.not.atomrigidcoordt)) then1060:          if (rigidinit.and.(.not.atomrigidcoordt)) then
1051:             xrigidcoords(1:degfreedoms) = x(1:degfreedoms)1061:             xrigidcoords(1:degfreedoms) = x(1:degfreedoms)
1052:             call transformrigidtoc(1, nrigidbody, x, xrigidcoords)1062:             call transformrigidtoc(1, nrigidbody, x, xrigidcoords)
1053:          endif1063:          endif
 1064:          !call ljpshift(x, grad, ereal, gradt, sect)
 1065:          !print *, 'ljpshift energy    : ', ereal
 1066:          !print *, 'ljpshift grad      : ', sum(grad(1:3*natoms))
 1067:          !do j1 = 1, 3*natoms
 1068:          !   print *, j1, grad(j1)
 1069:          !enddo!, grad(3)!sum(grad(1:3*natoms))
1054:          call ljpshift_new(x, ereal, grad, gradt)1070:          call ljpshift_new(x, ereal, grad, gradt)
 1071:          !print *, 'ljpshift_new energy: ', ereal
 1072:          !print *, 'ljpshift_new grad  : ', sum(grad(1:3*natoms))
 1073:          !do j1 = 1, 3*natoms
 1074:          !   print *, j1, grad(j1)
 1075:          !enddo!, grad(3)!sum(grad(1:3*natoms))
1055:          if (rigidinit.and.(.not.atomrigidcoordt)) then1076:          if (rigidinit.and.(.not.atomrigidcoordt)) then
1056:             x(degfreedoms+1:3*natoms)=0.0d01077:             x(degfreedoms+1:3*natoms)=0.0d0
1057:             x(1:degfreedoms)=xrigidcoords(1:degfreedoms)1078:             x(1:degfreedoms)=xrigidcoords(1:degfreedoms)
1058:             if (gradt) then1079:             if (gradt) then
1059:                call transformgrad(grad, xrigidcoords, xrigidgrad)1080:                call transformgrad(grad, xrigidcoords, xrigidgrad)
1060:                gradatoms(1:3*natoms) = grad(1:3*natoms)1081:                gradatoms(1:3*natoms) = grad(1:3*natoms)
1061:                grad(degfreedoms+1:3*natoms) = 0.0d01082:                grad(degfreedoms+1:3*natoms) = 0.0d0
1062:                grad(1:degfreedoms) = xrigidgrad(1:degfreedoms)1083:                grad(1:degfreedoms) = xrigidgrad(1:degfreedoms)
1063:             endif1084:             endif
1064:          endif1085:          endif
 1086:          !stop
1065:       ELSE1087:       ELSE
1066:          CALL LJPBIN(X, GRAD, EREAL, GRADT, SECT)1088:          CALL LJPBIN(X, GRAD, EREAL, GRADT, SECT)
1067:       END IF1089:       END IF
1068: 1090: 
1069:    ELSE IF (SOFT_SPHERE) THEN1091:    ELSE IF (SOFT_SPHERE) THEN
1070:       CALL SOFT_SPHERE_POT(X, GRAD, EREAL, GRADT, SECT)1092:       CALL SOFT_SPHERE_POT(X, GRAD, EREAL, GRADT, SECT)
1071: 1093: 
1072:    ELSE IF (LB2T) THEN1094:    ELSE IF (LB2T) THEN
1073:       CALL RADR(X, GRAD, EREAL, GRADT)1095:       CALL RADR(X, GRAD, EREAL, GRADT)
1074:       CALL LB2(X, GRAD, EREAL, GRADT)1096:       CALL LB2(X, GRAD, EREAL, GRADT)
1075: 1097: 
1076:    ELSE IF (LJCOULT) THEN1098:    ELSE IF (LJCOULT) THEN
1077:       CALL RADR(X, GRAD, EREAL, GRADT)1099:       CALL RADR(X, GRAD, EREAL, GRADT)
1078:       CALL LJCOUL(X, GRAD, EREAL, GRADT)1100:       CALL LJCOUL(X, GRAD, EREAL, GRADT)
1079: 1101: 
1080:   ELSE IF (LJ_GAUSST) THEN1102:   ELSE IF (LJ_GAUSST) THEN
 1103:       ! dj337: changed to test box derivatives
 1104:       print *, 'new pot call!!!!!'
 1105:       print *, 'frac coords: '
 1106:       do j = 1, natoms
 1107:          print *, j, x(3*j-2:3*j)
 1108:       enddo
1081:       CALL LJ_GAUSS(X, GRAD, EREAL, GRADT)1109:       CALL LJ_GAUSS(X, GRAD, EREAL, GRADT)
 1110:       print *, 'frac grad: '
 1111:       do j = 1, natoms
 1112:          print *, j, grad(3*j-2:3*j)
 1113:       enddo
 1114:       !print *, 'cart grad: '
 1115:       !do j = 1, natoms
 1116:       !   print *, j, grad2(3*j-2:3*j)
 1117:       !enddo
 1118:       !print *, 'val: ', val
 1119:       !print *, 'og grad: '
 1120:       !do j = 1, natoms
 1121:       !   print *, grad(3*j-2:3*j)
 1122:       !enddo
 1123:       !print *, 'og x: '
 1124:       !do j = 1, natoms-2
 1125:       !   print *, x(3*j-2:3*j)
 1126:       !enddo
 1127:       !print *, 'og grad: '
 1128:       !do j = 1, natoms
 1129:       !   print *, grad(3*j-2:3*j)
 1130:       !enddo
 1131:       !print *, 'og box params: '
 1132:       !print *, x(3*natoms-2:3*natoms)
 1133:       !print *, x(3*natoms-5:3*natoms-3)
 1134:       !print *, 'lj_gauss box derivs'
 1135:       !print *, grad(3*natoms-2:3*natoms)
 1136:       !print *, grad(3*natoms-5:3*natoms-3)
 1137:       !boxparams(4:6) = 1.57079633 !x(3*natoms-2:3*natoms)
 1138:       !boxparams(1:3) = x(3*natoms-5:3*natoms-3)
 1139:       !if (hithere) then
 1140:       !   print *, 'start: ', boxparamsgrad(1:6)
 1141:       !   call frac2cart(xfrac, x)
 1142:       !   call boxderiv_cart_tri(xfrac, grad2, xfracgrad)!, boxparams, boxparamsgrad)
 1143:       !   print *, 'cartesian coords:'
 1144:       !   do j = 1, natoms-2
 1145:       !      print *, j, xfrac(3*j-2:3*j)
 1146:       !   enddo
 1147:       !   print *, 'cartesian grad:'
 1148:       !   do j = 1, natoms-2
 1149:       !      print *, j, grad2(3*j-2:3*j)
 1150:       !   enddo
 1151:       !   !print *, 'final: ', boxparamsgrad(1:6)
 1152:       !   grad(1:3*natoms-6) = xfracgrad(1:3*natoms-6)
 1153:       !   grad(3*natoms-5:3*natoms) = 0.0d0
 1154:       !   grad(3*natoms-5:3*natoms-3) = boxparamsgrad(4:6)
 1155:       !   grad(3*natoms-2:3*natoms) = boxparamsgrad(1:3)
 1156:       !   print *, 'fake grad: '
 1157:       !   do j = 1, natoms-2
 1158:       !      print *, xfracgrad(3*j-2:3*j)
 1159:       !   enddo
 1160:       !   print *, boxparamsgrad(4:6)
 1161:       !   print *, boxparamsgrad(1:3)
 1162:       !   !print *, 'fake box params: '
 1163:       !   !print *, boxparams(1:3)
 1164:       !   !print *, boxparams(4:6)
 1165:       !   print *, 'rms: ', rms
 1166:       !   !print *, 'fake x: '
 1167:       !   !do j = 1, natoms
 1168:       !   !   print *, x(3*j-2:3*j)
 1169:       !   !enddo
 1170:       !endif
 1171:       stop
1082: 1172: 
1083:    ELSE IF (STOCKT) THEN1173:    ELSE IF (STOCKT) THEN
1084:       CALL RADR(X, GRAD, EREAL, GRADT)1174:       CALL RADR(X, GRAD, EREAL, GRADT)
1085:       CALL STOCK(X, GRAD, EREAL, GRADT)1175:       CALL STOCK(X, GRAD, EREAL, GRADT)
1086: 1176: 
1087: ! RBAA potentials1177: ! RBAA potentials
1088: 1178: 
1089:    ELSE IF (CAPBINT) THEN1179:    ELSE IF (CAPBINT) THEN
1090:       CALL CAPBIN (X, GRAD, EREAL, GRADT)1180:       CALL CAPBIN (X, GRAD, EREAL, GRADT)
1091: 1181: 
1384:          CALL LJCUT(X, GRAD, EREAL, GRADT, SECT)1474:          CALL LJCUT(X, GRAD, EREAL, GRADT, SECT)
1385:       ELSE1475:       ELSE
1386:          CALL LJ(X, GRAD, EREAL, GRADT, SECT)1476:          CALL LJ(X, GRAD, EREAL, GRADT, SECT)
1387:          IF (AXTELL) CALL AXT(NATOMS, X, GRAD, EREAL, GRADT, ZSTAR)1477:          IF (AXTELL) CALL AXT(NATOMS, X, GRAD, EREAL, GRADT, ZSTAR)
1388:       END IF1478:       END IF
1389:    END IF1479:    END IF
1390: !1480: !
1391: !  --------------- End of possible potentials - now add fields if required ------------------------------1481: !  --------------- End of possible potentials - now add fields if required ------------------------------
1392: !1482: !
1393: 1483: 
 1484:    !print *, 'end coords pot: '
 1485:    !do j1 = 1, 2*nrigidbody
 1486:    !   print *, x(3*j1-2:3*j1)
 1487:    !enddo
 1488: 
1394: ! dj337: compute box derivatives1489: ! dj337: compute box derivatives
1395:     if (boxderivt) then1490:     if (boxderivt) then
1396:        ! convert from absolute to fractional coordinates1491:        ! convert from absolute to fractional coordinates
1397:        if (.not.(rigidinit.and.atomrigidcoordt)) then1492:        if (.not.(rigidinit.and.atomrigidcoordt)) then
1398:           call cart2frac(x, xfrac)1493:           call cart2frac(x, xfrac)
1399:           x(1:3*natoms) = xfrac(1:3*natoms)1494:           x(1:3*natoms) = xfrac(1:3*natoms)
1400:        endif1495:        endif
1401:        ! convert gradient from absolute to fractional1496:        ! convert gradient from absolute to fractional
1402:        if (gradt) then1497:        if (gradt) then
1403:           call boxderiv(grad, gradfrac)1498:           call boxderiv(grad, gradfrac)
1404:           grad(1:3*natoms) = gradfrac(1:3*natoms)1499:           grad(1:3*natoms) = gradfrac(1:3*natoms)
1405:        endif1500:        endif
1406:     endif1501:     endif
1407: 1502: 
 1503:    !print *, 'end2 coords pot: '
 1504:    !do j1 = 1, 2*nrigidbody
 1505:    !   print *, x(3*j1-2:3*j1)
 1506:    !enddo
 1507: 
1408:    IF (RIGIDCONTOURT) THEN1508:    IF (RIGIDCONTOURT) THEN
1409: !       sf344> sanity check, since this keyword works only for a pair of rigid particles1509: !       sf344> sanity check, since this keyword works only for a pair of rigid particles
1410:       IF (NATOMS.NE.4) THEN1510:       IF (NATOMS.NE.4) THEN
1411:          WRITE(MYUNIT,'(A)') 'potential> ERROR - RIGIDCONTOUR works only for a PAIR of rigid bodies, exiting'1511:          WRITE(MYUNIT,'(A)') 'potential> ERROR - RIGIDCONTOUR works only for a PAIR of rigid bodies, exiting'
1412:          STOP1512:          STOP
1413:       ELSE1513:       ELSE
1414:          CALL RIGIDCONTOUR()1514:          CALL RIGIDCONTOUR()
1415:          WRITE(MYUNIT,'(A)') 'potential> finished generating the potential energy matrix, exiting GMIN'1515:          WRITE(MYUNIT,'(A)') 'potential> finished generating the potential energy matrix, exiting GMIN'
1416:          STOP1516:          STOP
1417:       END IF1517:       END IF


r33132/saveit.f 2017-08-07 15:31:05.053716358 +0100 r33131/saveit.f 2017-08-07 15:31:07.717751796 +0100
 20:       USE COMMONS 20:       USE COMMONS
 21:       use QMODULE 21:       use QMODULE
 22:       IMPLICIT NONE 22:       IMPLICIT NONE
 23:  23: 
 24:  24: 
 25:       INTEGER J1, J2, J3, NP, NQTOT, CSMIT  25:       INTEGER J1, J2, J3, NP, NQTOT, CSMIT 
 26:       DOUBLE PRECISION EREAL,P(3*NATOMS), AVVAL, CSMRMS 26:       DOUBLE PRECISION EREAL,P(3*NATOMS), AVVAL, CSMRMS
 27:  27: 
 28:       COMMON /TOT/ NQTOT 28:       COMMON /TOT/ NQTOT
 29:      29:     
  30:       !print *, 'saveit: ', qminp(1, 1:5)
  31:       !print *, 'saveitbp: ', boxq(1, 1:6)
 30: ! 32: !
 31: ! Save the lowest minimum for each size with GCBH. 33: ! Save the lowest minimum for each size with GCBH.
 32: ! 34: !
 33:       IF (GCBHT) THEN 35:       IF (GCBHT) THEN
 34: !        IF (EREAL.LT.QENERGIES(NATOMS)) THEN 36: !        IF (EREAL.LT.QENERGIES(NATOMS)) THEN
 35:          IF (FEBH_POT_ENE+NATOMS*GCMU.LT.QPE(NATOMS)) THEN 37:          IF (FEBH_POT_ENE+NATOMS*GCMU.LT.QPE(NATOMS)) THEN
 36:             QENERGIES(NATOMS)=EREAL 38:             QENERGIES(NATOMS)=EREAL
 37:             QPE(NATOMS)=FEBH_POT_ENE+NATOMS*GCMU 39:             QPE(NATOMS)=FEBH_POT_ENE+NATOMS*GCMU
 38:             DO J2=1,3*NATOMS 40:             DO J2=1,3*NATOMS
 39:                QCOORDINATES(NATOMS,J2)=P(J2) 41:                QCOORDINATES(NATOMS,J2)=P(J2)
 40:             ENDDO 42:             ENDDO
 41:          ENDIF 43:          ENDIF
 42:       ENDIF 44:       ENDIF
 43: C 45: C
 44: C  Save the lowest NSAVE distinguishable configurations. 46: C  Save the lowest NSAVE distinguishable configurations.
 45: C 47: C
 46: !      WRITE(MYUNIT,'(A,12E15.7)') 'EREAL,ECONV,QMIN',EREAL,ECONV,(QMIN(J1),J1=1,NSAVE) 48: !      WRITE(MYUNIT,'(A,12E15.7)') 'EREAL,ECONV,QMIN',EREAL,ECONV,(QMIN(J1),J1=1,NSAVE)
 47:       DO J1=1,NSAVE 49:       DO J1=1,NSAVE
  50:          !print *, 'saveit 48'
 48:          IF (DABS(EREAL-QMIN(J1)).LT.ECONV) THEN 51:          IF (DABS(EREAL-QMIN(J1)).LT.ECONV) THEN
 49: C 52: C
 50: C  These are probably the same - but just to make sure we save the  53: C  These are probably the same - but just to make sure we save the 
 51: C  lowest. 54: C  lowest.
 52: C 55: C
 53:             IF (EREAL.LT.QMIN(J1)) THEN 56:             IF (EREAL.LT.QMIN(J1)) THEN
  57:                !print *, 'saveit 55'
 54:                QMINNATOMS(J1)=NATOMS 58:                QMINNATOMS(J1)=NATOMS
 55:                QMIN(J1)=EREAL 59:                QMIN(J1)=EREAL
 56:                DO J2=1,3*NATOMS 60:                DO J2=1,3*NATOMS
 57:                   QMINP(J1,J2)=P(J2) 61:                   QMINP(J1,J2)=P(J2)
 58:                ENDDO 62:                ENDDO
 59:                ! dj337: save box parameters 63:                !print *, 'gsaving boxparams'
 60:                if (boxderivt) boxq(j1,1:6) = box_params(1:6) 64:                if (boxderivt) boxq(j1,1:6) = box_params(1:6)
 61:             ENDIF 65:             ENDIF
 62:             GOTO 10 66:             GOTO 10
 63:          ENDIF 67:          ENDIF
 64:          IF (EREAL.LT.QMIN(J1)) THEN 68:          IF (EREAL.LT.QMIN(J1)) THEN
  69:             !print *, 'saveit 67'
 65:  70: 
 66:             J2=NSAVE 71:             J2=NSAVE
 67: 20          CONTINUE 72: 20          CONTINUE
  73:             !print *, 'saveit 71' 
 68:             IF (NSAVE.GT.1) THEN 74:             IF (NSAVE.GT.1) THEN
  75:                !print *, 'saveit 73'
 69:                QMIN(J2)=QMIN(J2-1) 76:                QMIN(J2)=QMIN(J2-1)
 70:                QMINNATOMS(J2)=QMINNATOMS(J2-1) 77:                QMINNATOMS(J2)=QMINNATOMS(J2-1)
 71:                FF(J2)=FF(J2-1) 78:                FF(J2)=FF(J2-1)
 72:                NPCALL_QMIN(J2)=NPCALL_QMIN(J2-1) 79:                NPCALL_QMIN(J2)=NPCALL_QMIN(J2-1)
 73:                DO J3=1,3*QMINNATOMS(J2-1) 80:                DO J3=1,3*QMINNATOMS(J2-1)
 74:                   QMINP(J2,J3)=QMINP(J2-1,J3) 81:                   QMINP(J2,J3)=QMINP(J2-1,J3)
 75:                ENDDO 82:                ENDDO
 76:                ! dj337: update saved box parameters 
 77:                if (boxderivt) boxq(j2,1:6) = boxq(j2-1,1:6) 83:                if (boxderivt) boxq(j2,1:6) = boxq(j2-1,1:6)
 78:  84: 
 79:                J2=J2-1 85:                J2=J2-1
 80:                IF (J2.GE.J1+1) GOTO 20 86:                IF (J2.GE.J1+1) GOTO 20
 81:             ENDIF 87:             ENDIF
  88:             !print *, 'saveit2: ', qminp(1, 1:5)
  89:             !print *, 'saveitbp2: ', boxq(1, 1:6)
 82:             QMIN(J1)=EREAL 90:             QMIN(J1)=EREAL
 83:             QMINNATOMS(J1)=NATOMS 91:             QMINNATOMS(J1)=NATOMS
 84:             FF(J1)=NQ(NP) 92:             FF(J1)=NQ(NP)
 85:             NPCALL_QMIN(J1)=NPCALL 93:             NPCALL_QMIN(J1)=NPCALL
 86:             DO J2=1,3*QMINNATOMS(J1) 94:             DO J2=1,3*QMINNATOMS(J1)
 87:                QMINP(J1,J2)=P(J2) 95:                QMINP(J1,J2)=P(J2)
 88:             ENDDO 96:             ENDDO
 89:             ! dj337: save box parameters 
 90:             if (boxderivt) boxq(j1,1:6) = box_params(1:6) 97:             if (boxderivt) boxq(j1,1:6) = box_params(1:6)
  98:             !print *, 'saveit3: ', qminp(1, 1:5)
  99:             !print *, 'saveitbp3: ', boxq(1, 1:6)
 91: 100: 
 92:             GOTO 10101:             GOTO 10
 93:          ENDIF102:          ENDIF
 94:       ENDDO103:       ENDDO
 95: 104: 
 96: 10    CONTINUE105: 10    CONTINUE
 97: 106: 
 98: 107: 
 99: ! beh35 > adding stuff to accommodate for mlowest108: ! beh35 > adding stuff to accommodate for mlowest
100: 109: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0