hdiff output

r32456/benzgenrigid.f90 2017-05-04 10:30:15.730752365 +0100 r32455/benzgenrigid.f90 2017-05-04 10:30:16.222758732 +0100
233:                                 + 1.5D0*COSTB*COSTB - 1.D0)233:                                 + 1.5D0*COSTB*COSTB - 1.D0)
234:                         ! ENERGY1 is energy due to short-range anisotropic interactions234:                         ! ENERGY1 is energy due to short-range anisotropic interactions
235:                         ! calculate vertical shift for first term235:                         ! calculate vertical shift for first term
236:                         EXPFCT  = KKJ*DEXP(-ALPHACC*(ABSRIJ - RHOCC))236:                         EXPFCT  = KKJ*DEXP(-ALPHACC*(ABSRIJ - RHOCC))
237:                         VSHIFT1 = KKJ*DEXP(-ALPHACC*(EWALDREALC - RHOCC))237:                         VSHIFT1 = KKJ*DEXP(-ALPHACC*(EWALDREALC - RHOCC))
238:                         ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1238:                         ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
239:                         ! ENERGY2 is energy due to damped dispersion239:                         ! ENERGY2 is energy due to damped dispersion
240:                         ! calculate vertical shift for second term240:                         ! calculate vertical shift for second term
241:                         VSHIFT2 = DC6CC*DMPFCT_SHIFT/(EWALDREALC**6)241:                         VSHIFT2 = DC6CC*DMPFCT_SHIFT/(EWALDREALC**6)
242:                         !print *, 'energy: ', dc6cc*dmpfct*r6242:                         !print *, 'energy: ', dc6cc*dmpfct*r6
243:                         ENERGY2 = ENERGY2 - DC6CC*DMPFCT*R6 + VSHIFT2243:                         ENERGY2 = ENERGY2 - DC6CC*DMPFCT*R6 - VSHIFT2
244:    244:    
245:                         IF (GTEST) THEN245:                         IF (GTEST) THEN
246:    246:    
247:                            ! DVDR is derivative of dispersion damping factor energy with factor 1/Rab247:                            ! DVDR is derivative of dispersion damping factor energy with factor 1/Rab
248:                            DVDR    = 6.D0*DC6CC*R6*R2*DMPFCT - DC6CC*R6*DDMPDR 248:                            DVDR    = 6.D0*DC6CC*R6*R2*DMPFCT - DC6CC*R6*DDMPDR 
249:                            !print *, 'grad: ', dvdr249:                            !print *, 'grad: ', dvdr
250:                            ! FRIJ is derivative of ENERGY1 wrt r_ij with factor 1/Rab250:                            ! FRIJ is derivative of ENERGY1 wrt r_ij with factor 1/Rab
251:                            FRIJ(:) = ALPHACC*EXPFCT*(-NR(:) + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADR(:) &251:                            FRIJ(:) = ALPHACC*EXPFCT*(-NR(:) + (RHOCC10 + 3.D0*RHOCC20*COSTA)*DCADR(:) &
252:                                    + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDR(:))252:                                    + (RHOCC10 + 3.D0*RHOCC20*COSTB)*DCBDR(:))
253:                            ! TIJ is derivative of ENERGY1 wrt pi with factor 1/Rab253:                            ! TIJ is derivative of ENERGY1 wrt pi with factor 1/Rab
262:                      ! calculate if I and J are both hydorgens262:                      ! calculate if I and J are both hydorgens
263:                      ELSEIF (I > NCARBON .AND. J > NCARBON) THEN263:                      ELSEIF (I > NCARBON .AND. J > NCARBON) THEN
264:    264:    
265:                         RHOHH  = RHOHH0 + RHOHH10*(COSTA + COSTB) + RHOHH20*(1.5D0*COSTA*COSTA      &265:                         RHOHH  = RHOHH0 + RHOHH10*(COSTA + COSTB) + RHOHH20*(1.5D0*COSTA*COSTA      &
266:                                + 1.5D0*COSTB*COSTB - 1.D0) 266:                                + 1.5D0*COSTB*COSTB - 1.D0) 
267:                         EXPFCT  = KKJ*DEXP(-ALPHAHH*(ABSRIJ - RHOHH))267:                         EXPFCT  = KKJ*DEXP(-ALPHAHH*(ABSRIJ - RHOHH))
268:                         VSHIFT1 = KKJ*DEXP(-ALPHAHH*(EWALDREALC - RHOHH))268:                         VSHIFT1 = KKJ*DEXP(-ALPHAHH*(EWALDREALC - RHOHH))
269:                         ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1269:                         ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
270:                         VSHIFT2 = DC6HH*DMPFCT_SHIFT/(EWALDREALC**6)270:                         VSHIFT2 = DC6HH*DMPFCT_SHIFT/(EWALDREALC**6)
271:                         !print *, 'energy: ', dc6hh*dmpfct*r6271:                         !print *, 'energy: ', dc6hh*dmpfct*r6
272:                         ENERGY2 = ENERGY2 - DC6HH*DMPFCT*R6 + VSHIFT2272:                         ENERGY2 = ENERGY2 - DC6HH*DMPFCT*R6 - VSHIFT2
273:    273:    
274:                         IF (GTEST) THEN274:                         IF (GTEST) THEN
275:    275:    
276:                            DVDR    = 6.D0*DC6HH*R6*R2*DMPFCT - DC6HH*R6*DDMPDR 276:                            DVDR    = 6.D0*DC6HH*R6*R2*DMPFCT - DC6HH*R6*DDMPDR 
277:                            !print *, 'grad: ', dvdr277:                            !print *, 'grad: ', dvdr
278:                            FRIJ(:) = ALPHAHH*EXPFCT*(-NR(:) + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADR(:) &278:                            FRIJ(:) = ALPHAHH*EXPFCT*(-NR(:) + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADR(:) &
279:                                    + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDR(:))279:                                    + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDR(:))
280:                            TIJ(:)  = ALPHAHH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPI(:) &280:                            TIJ(:)  = ALPHAHH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPI(:) &
281:                                    + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPI(:))281:                                    + (RHOHH10 + 3.D0*RHOHH20*COSTB)*DCBDPI(:))
282:                            TJI(:)  = ALPHAHH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPJ(:) &282:                            TJI(:)  = ALPHAHH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOHH10 + 3.D0*RHOHH20*COSTA)*DCADPJ(:) &
287:                      ! calculate if I is carbon and J is hydrogen287:                      ! calculate if I is carbon and J is hydrogen
288:                      ELSE IF (I <= NCARBON .AND. J > NCARBON) THEN 288:                      ELSE IF (I <= NCARBON .AND. J > NCARBON) THEN 
289:    289:    
290:                         RHOCH  = RHOCH0 + RHOC10H*COSTA + RHOCH10*COSTB + RHOC20H*(1.5D0*COSTA*COSTA &290:                         RHOCH  = RHOCH0 + RHOC10H*COSTA + RHOCH10*COSTB + RHOC20H*(1.5D0*COSTA*COSTA &
291:                                - 0.5D0) + RHOCH20*(1.5D0*COSTB*COSTB - 0.5D0)291:                                - 0.5D0) + RHOCH20*(1.5D0*COSTB*COSTB - 0.5D0)
292:                         EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))292:                         EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))
293:                         VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))293:                         VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))
294:                         ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1294:                         ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
295:                         VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)295:                         VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)
296:                         !print *, 'energy: ', dc6ch*dmpfct*r6296:                         !print *, 'energy: ', dc6ch*dmpfct*r6
297:                         ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2297:                         ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 - VSHIFT2
298:    298:    
299:                         IF (GTEST) THEN299:                         IF (GTEST) THEN
300:                      300:                      
301:                            DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 301:                            DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 
302:                            !print *, 'grad: ', dvdr302:                            !print *, 'grad: ', dvdr
303:                            FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADR(:) &303:                            FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADR(:) &
304:                                    + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDR(:))304:                                    + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDR(:))
305:                            TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPI(:) &305:                            TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPI(:) &
306:                                    + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPI(:))306:                                    + (RHOCH10 + 3.D0*RHOCH20*COSTB)*DCBDPI(:))
307:                            TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPJ(:) &307:                            TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOC10H + 3.D0*RHOC20H*COSTA)*DCADPJ(:) &
311:    311:    
312:                      ELSE !IF(I > NCARBON .AND. J <= NCARBON) THEN312:                      ELSE !IF(I > NCARBON .AND. J <= NCARBON) THEN
313:    313:    
314:                         RHOCH  = RHOCH0 + RHOCH10*COSTA + RHOC10H*COSTB + RHOCH20*(1.5D0*COSTA*COSTA &314:                         RHOCH  = RHOCH0 + RHOCH10*COSTA + RHOC10H*COSTB + RHOCH20*(1.5D0*COSTA*COSTA &
315:                                - 0.5D0) + RHOC20H*(1.5D0*COSTB*COSTB - 0.5D0)315:                                - 0.5D0) + RHOC20H*(1.5D0*COSTB*COSTB - 0.5D0)
316:                         EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))316:                         EXPFCT  = KKJ*DEXP(-ALPHACH*(ABSRIJ - RHOCH))
317:                         VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))317:                         VSHIFT1 = KKJ*DEXP(-ALPHACH*(EWALDREALC - RHOCH))
318:                         ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1318:                         ENERGY1 = ENERGY1 + EXPFCT - VSHIFT1
319:                         VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)319:                         VSHIFT2 = DC6CH*DMPFCT_SHIFT/(EWALDREALC**6)
320:                         !print *, 'energy: ', dc6ch*dmpfct*r6320:                         !print *, 'energy: ', dc6ch*dmpfct*r6
321:                         ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 + VSHIFT2321:                         ENERGY2 = ENERGY2 - DC6CH*DMPFCT*R6 - VSHIFT2
322:    322:    
323:                         IF (GTEST) THEN323:                         IF (GTEST) THEN
324:    324:    
325:                            DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 325:                            DVDR    = 6.D0*DC6CH*R6*R2*DMPFCT - DC6CH*R6*DDMPDR 
326:                            !print *, 'grad: ', dvdr326:                            !print *, 'grad: ', dvdr
327:                            FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADR(:) &327:                            FRIJ(:) = ALPHACH*EXPFCT*(-NR(:) + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADR(:) &
328:                                    + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDR(:))328:                                    + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDR(:))
329:                            TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPI(:) &329:                            TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPI(:) &
330:                                    + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPI(:))330:                                    + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPI(:))
331:                            TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPJ(:) &331:                            TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPJ(:) &


r32456/potential.f90 2017-05-04 10:30:15.970755469 +0100 r32455/potential.f90 2017-05-04 10:30:16.458761787 +0100
 39:    USE MBPOLMOD, ONLY: MBPOL 39:    USE MBPOLMOD, ONLY: MBPOL
 40:    USE SWMOD, ONLY: SWTYPE 40:    USE SWMOD, ONLY: SWTYPE
 41:    USE MODCUDALBFGS, ONLY: CUDA_ENEGRAD_WRAPPER, CUDA_NUMERICAL_HESS 41:    USE MODCUDALBFGS, ONLY: CUDA_ENEGRAD_WRAPPER, CUDA_NUMERICAL_HESS
 42:    USE GAUSS_MOD, ONLY: GFIELD 42:    USE GAUSS_MOD, ONLY: GFIELD
 43:    USE RAD_MOD, ONLY: RAD, RADR 43:    USE RAD_MOD, ONLY: RAD, RADR
 44:    USE PREC, ONLY: INT32, REAL64 44:    USE PREC, ONLY: INT32, REAL64
 45:    USE TWIST_MOD, ONLY: TWISTT, TWIST 45:    USE TWIST_MOD, ONLY: TWISTT, TWIST
 46:    USE CONVEX_POLYHEDRA_MODULE, ONLY: CONVEX_POLYHEDRA 46:    USE CONVEX_POLYHEDRA_MODULE, ONLY: CONVEX_POLYHEDRA
 47:    USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS 47:    USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS
 48:    USE EWALD 48:    USE EWALD
 49:    USE BOXDERIV 
 50:    IMPLICIT NONE 49:    IMPLICIT NONE
 51: ! Arguments 50: ! Arguments
 52: ! TODO: Work out intents 51: ! TODO: Work out intents
 53: ! TODO: Fix array dimensions? 52: ! TODO: Fix array dimensions?
 54:    REAL(REAL64)               :: X(*) 53:    REAL(REAL64)               :: X(*)
 55:    REAL(REAL64)               :: GRAD(*) 54:    REAL(REAL64)               :: GRAD(*)
 56:    REAL(REAL64)               :: EREAL 55:    REAL(REAL64)               :: EREAL
 57:    LOGICAL, INTENT(IN)        :: GRADT 56:    LOGICAL, INTENT(IN)        :: GRADT
 58:    LOGICAL, INTENT(IN)        :: SECT 57:    LOGICAL, INTENT(IN)        :: SECT
 59: ! Variables 58: ! Variables
 60: ! hk286 > generalised rigid body additions 59: ! hk286 > generalised rigid body additions
 61:    REAL(REAL64)               :: XCOORDS(3*NATOMS), GRADATOMS(3*NATOMS), & 60:    REAL(REAL64)               :: XCOORDS(3*NATOMS), GRADATOMS(3*NATOMS), &
 62:                                  XRIGIDCOORDS(DEGFREEDOMS), XRIGIDGRAD(DEGFREEDOMS), XFRAC(3*NATOMS) 61:                                  XRIGIDCOORDS(DEGFREEDOMS), XRIGIDGRAD(DEGFREEDOMS)
 63:    REAL(REAL64), ALLOCATABLE  :: TEMPHESS(:) 62:    REAL(REAL64), ALLOCATABLE  :: TEMPHESS(:)
 64:    REAL(REAL64)               :: GRAD1(3*NATOMS) 63:    REAL(REAL64)               :: GRAD1(3*NATOMS)
 65:    INTEGER(INT32)             :: I, J, K 64:    INTEGER(INT32)             :: I, J, K
 66:    REAL(REAL64)               :: XRIGIDHESS(DEGFREEDOMS, DEGFREEDOMS) 65:    REAL(REAL64)               :: XRIGIDHESS(DEGFREEDOMS, DEGFREEDOMS)
 67:    LOGICAL                    :: FTEST, EVAP, COMPON, YESNO, GUIDET, EVAPREJECT 66:    LOGICAL                    :: FTEST, EVAP, COMPON, YESNO, GUIDET, EVAPREJECT
 68:    INTEGER(INT32)             :: J1, J2, J3, CSMIT 67:    INTEGER(INT32)             :: J1, J2, J3, CSMIT
 69:    CHARACTER                  :: FNAME*80, DUMM*4 68:    CHARACTER                  :: FNAME*80, DUMM*4
 70:    REAL(REAL64)               :: DUMMY2, GEMAX, RMAT(3, 3), XTEMP(3*NATOMS), AVVAL, & 69:    REAL(REAL64)               :: DUMMY2, GEMAX, RMAT(3, 3), XTEMP(3*NATOMS), AVVAL, &
 71:                                  DUMMY(3*NATOMS), DIST2, QE, QX, AA(3), SAVECSMNORM, & 70:                                  DUMMY(3*NATOMS), DIST2, QE, QX, AA(3), SAVECSMNORM, &
 72:                                  CSMRMS, CSMGRAD(3), SAVECSMPMAT(3, 3), & 71:                                  CSMRMS, CSMGRAD(3), SAVECSMPMAT(3, 3), &
 73:                                  SAVECSMIMAGES(3*NATOMS*CSMGPINDEX), BOXPARAMS(6) 72:                                  SAVECSMIMAGES(3*NATOMS*CSMGPINDEX)
 74:    CHARACTER(LEN=3)           :: CSMGPSAVE 73:    CHARACTER(LEN=3)           :: CSMGPSAVE
 75:    INTEGER(INT32)             :: CSMGPINDEXSAVE 74:    INTEGER(INT32)             :: CSMGPINDEXSAVE
 76:    REAL(REAL64)               :: PTGPSAVE(3, 3, 2*CSMGPINDEX), CSMNORMSAVE, ENERGY, VNEW(3*NATOMS) 75:    REAL(REAL64)               :: PTGPSAVE(3, 3, 2*CSMGPINDEX), CSMNORMSAVE, ENERGY, VNEW(3*NATOMS)
 77:    INTEGER(INT32)             :: BRUN 76:    INTEGER(INT32)             :: BRUN
 78:    INTEGER(INT32)             :: NQTOT 77:    INTEGER(INT32)             :: NQTOT
 79:    LOGICAL                    :: SOCOUPLE, GUIDECHANGET, CSMDOGUIDET, COMTEST 78:    LOGICAL                    :: SOCOUPLE, GUIDECHANGET, CSMDOGUIDET, COMTEST
 80: ! khs26 > AMBER12 energy decomposition, type defined in AMBER12 interface 79: ! khs26 > AMBER12 energy decomposition, type defined in AMBER12 interface
 81:    TYPE(POT_ENE_REC_C)        :: AMBER12_ENERGY_DECOMP 80:    TYPE(POT_ENE_REC_C)        :: AMBER12_ENERGY_DECOMP
 82: !Used only in commented sections 81: !Used only in commented sections
 83:    REAL(REAL64)               :: XG, YG, ZG, GRADLJ(3*NATOMS), EREALLJ, GRADMF(3*NATOMS), EREALMF, TERMLJ, TERMMF, & 82:    REAL(REAL64)               :: XG, YG, ZG, GRADLJ(3*NATOMS), EREALLJ, GRADMF(3*NATOMS), EREALMF, TERMLJ, TERMMF, &
324: !         WRITE(*,'(I5, 4F20.10)') J1, GRAD(J1), EPLUS, EMINUS, ABS(EPLUS-EMINUS)323: !         WRITE(*,'(I5, 4F20.10)') J1, GRAD(J1), EPLUS, EMINUS, ABS(EPLUS-EMINUS)
325: !     END DO324: !     END DO
326: !325: !
327: 326: 
328: ! dj337: pahagenrigid for benzene with Ewald327: ! dj337: pahagenrigid for benzene with Ewald
329:    ELSE IF (BENZRIGIDEWALDT) THEN328:    ELSE IF (BENZRIGIDEWALDT) THEN
330:       CALL BENZGENRIGIDEWALD(X, GRAD, EREAL, GRADT)329:       CALL BENZGENRIGIDEWALD(X, GRAD, EREAL, GRADT)
331:       if (rigidinit .and. (.not.atomrigidcoordt)) then330:       if (rigidinit .and. (.not.atomrigidcoordt)) then
332:          xrigidcoords(1:degfreedoms) = x(1:degfreedoms)331:          xrigidcoords(1:degfreedoms) = x(1:degfreedoms)
333:          xrigidgrad(1:degfreedoms) = grad(1:degfreedoms)332:          xrigidgrad(1:degfreedoms) = grad(1:degfreedoms)
334:          call transformrigidtoc(1, nrigidbody, xcoords, xrigidcoords) 
335:       endif333:       endif
336: ! check analytical and numerical gradients334: ! check analytical and numerical gradients
337: !      DIFF=1.0D-10335: !      DIFF=1.0D-10
338: !      WRITE(*, *) 'analytic and numerical gradients:'336: !      WRITE(*, *) 'analytic and numerical gradients:'
339: !      print *, 'diff: ', diff337: !      print *, 'diff: ', diff
340: !      IF (ATOMRIGIDCOORDT) THEN ! if input is cartesian338: !      IF (ATOMRIGIDCOORDT) THEN ! if input is cartesian
341: !         XRIGIDCOORDS(:) = 0.D0339: !         XRIGIDCOORDS(:) = 0.D0
342: !         CALL TRANSFORMCTORIGID(X, XRIGIDCOORDS)340: !         CALL TRANSFORMCTORIGID(X, XRIGIDCOORDS)
343: !         X(1:DEGFREEDOMS) = XRIGIDCOORDS(1:DEGFREEDOMS)341: !         X(1:DEGFREEDOMS) = XRIGIDCOORDS(1:DEGFREEDOMS)
344: !      ENDIF342: !      ENDIF
1018:       CALL RAD(X, GRAD, EREAL, GRADT)1016:       CALL RAD(X, GRAD, EREAL, GRADT)
1019:       CALL GLJ(X, GRAD, EREAL, GRADT)1017:       CALL GLJ(X, GRAD, EREAL, GRADT)
1020: 1018: 
1021:    ELSE IF (BINARY) THEN1019:    ELSE IF (BINARY) THEN
1022:       IF (SHIFTCUT) THEN1020:       IF (SHIFTCUT) THEN
1023:          if (rigidinit.and.(.not.atomrigidcoordt)) then1021:          if (rigidinit.and.(.not.atomrigidcoordt)) then
1024:             xrigidcoords(1:degfreedoms) = x(1:degfreedoms)1022:             xrigidcoords(1:degfreedoms) = x(1:degfreedoms)
1025:             call transformrigidtoc(1, nrigidbody, x, xrigidcoords)1023:             call transformrigidtoc(1, nrigidbody, x, xrigidcoords)
1026:          endif1024:          endif
1027:          call ljpshift(x, grad, ereal, gradt, sect)1025:          call ljpshift(x, grad, ereal, gradt, sect)
1028:          boxparams = (/2.0D0, 1.0D0, 2.0D0, 3.0d0, 1.0d0, 2.0d0/) 
1029:          call cart2frac(x, xfrac, boxparams) 
1030:          print *, 'coords: ', x(:3*natoms) 
1031:          print *, 'frac: ', xfrac(:3*natoms) 
1032:          if (rigidinit.and.(.not.atomrigidcoordt)) then1026:          if (rigidinit.and.(.not.atomrigidcoordt)) then
1033:             x(degfreedoms+1:3*natoms)=0.0d01027:             x(degfreedoms+1:3*natoms)=0.0d0
1034:             x(1:degfreedoms)=xrigidcoords(1:degfreedoms)1028:             x(1:degfreedoms)=xrigidcoords(1:degfreedoms)
1035:             if (gradt) then1029:             if (gradt) then
1036:                call transformgrad(grad, xrigidcoords, xrigidgrad)1030:                call transformgrad(grad, xrigidcoords, xrigidgrad)
1037:                gradatoms(1:3*natoms) = grad(1:3*natoms)1031:                gradatoms(1:3*natoms) = grad(1:3*natoms)
1038:                grad(degfreedoms+1:3*natoms) = 0.0d01032:                grad(degfreedoms+1:3*natoms) = 0.0d0
1039:                grad(1:degfreedoms) = xrigidgrad(1:degfreedoms)1033:                grad(1:degfreedoms) = xrigidgrad(1:degfreedoms)
1040:             endif1034:             endif
1041:          endif1035:          endif
1042:          stop 
1043:       ELSE1036:       ELSE
1044:          CALL LJPBIN(X, GRAD, EREAL, GRADT, SECT)1037:          CALL LJPBIN(X, GRAD, EREAL, GRADT, SECT)
1045:       END IF1038:       END IF
1046: 1039: 
1047:    ELSE IF (SOFT_SPHERE) THEN1040:    ELSE IF (SOFT_SPHERE) THEN
1048:       CALL SOFT_SPHERE_POT(X, GRAD, EREAL, GRADT, SECT)1041:       CALL SOFT_SPHERE_POT(X, GRAD, EREAL, GRADT, SECT)
1049: 1042: 
1050:    ELSE IF (LB2T) THEN1043:    ELSE IF (LB2T) THEN
1051:       CALL RADR(X, GRAD, EREAL, GRADT)1044:       CALL RADR(X, GRAD, EREAL, GRADT)
1052:       CALL LB2(X, GRAD, EREAL, GRADT)1045:       CALL LB2(X, GRAD, EREAL, GRADT)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0