hdiff output

r33133/benzgenrigid.f90 2017-08-07 15:30:35.989329718 +0100 r33132/benzgenrigid.f90 2017-08-07 15:30:38.193359035 +0100
  1: ! -----------------------------------------------------------------------------  1: ! -----------------------------------------------------------------------------
  2: ! dj337: Anisotropic potential for periodic benzene systems.  2: ! dj337
  3: ! Potential from:  3: ! Anisotropic potential for periodic benzene systems. Potential from:
  4: ! Totton, TS, Misquitta, AJ, Kraft, M; J. Chem. Theory Comput., 2010, 6, 683-695.  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.
  5:   7: 
  6: ! Long-range electrostatic interactions computed using Ewald summation. 
  7: ! Implemented within the GENRIGID framework for rigid bodies.  8: ! Implemented within the GENRIGID framework for rigid bodies.
  8:   9: 
  9: ! If BOXDERIV keyword is used, energy gradients with respect to cell parameters 10: ! Gradients wrt lattice parameters are computed when BOXDERIV keyword used. This
 10: ! are computed to allow for unit cell optimization during the simulation. 11: ! subroutine is for triclinic systems. Use BENZGENRIGIDEWALD_ORTHO for orthorhombic
 11: ! This subroutine is for triclinic cell systems, use the subroutine  12: ! cells.
 12: ! BENZGENRIGIDEWALD_ORTHO for orthorhombic cell systems. 
 13: ! ----------------------------------------------------------------------------- 13: ! -----------------------------------------------------------------------------
 14:  14: 
 15:       SUBROUTINE BENZGENRIGIDEWALD(X, G, ENERGY, GTEST) 15:       SUBROUTINE BENZGENRIGIDEWALD(X, G, ENERGY, GTEST)
 16:  16: 
 17:       USE COMMONS, ONLY: NATOMS, NCARBON, RBSTLA, RHOCC0, RHOCC10, RHOCC20, & 17:       USE COMMONS, ONLY: NATOMS, NCARBON, RBSTLA, RHOCC0, RHOCC10, RHOCC20, &
 18:      &                   RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, RHOCH10, RHOC20H, & 18:      &                   RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, RHOCH10, RHOC20H, &
 19:      &                   RHOCH20, ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, & 19:      &                   RHOCH20, ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, &
 20:      &                   EWALDREALC, BOX_PARAMS, BOX_PARAMSGRAD 20:      &                   EWALDREALC, BOX_PARAMS, BOX_PARAMSGRAD
 21:  21: 
 22:       ! adapted to the genrigid framework 22:       ! adapted to the genrigid framework
 47:       DOUBLE PRECISION :: RHOCC, RHOHH, RHOCH, COSTA, COSTB, DMPFCT, DDMPDR, EXPFCT  47:       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) 48:       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) 49:       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) 50:       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) 51:       double precision :: reciplatvec(3,3), reciplatvec_grad(3,3,6)
 52:       DOUBLE PRECISION, PARAMETER :: B = 1.6485D0 52:       DOUBLE PRECISION, PARAMETER :: B = 1.6485D0
 53:       double precision, parameter :: pi = 3.141592654d0 53:       double precision, parameter :: pi = 3.141592654d0
 54:       integer, parameter          :: image_cutoff = 5 54:       integer, parameter          :: image_cutoff = 5
 55:       LOGICAL          :: GTEST, keep_angles 55:       LOGICAL          :: GTEST, keep_angles
 56:  56: 
 57:       ! check if combination of tricinlic cell angles is realistic 57:       ! check if combination of triclinic cell angles is realistic
 58:       if (boxderivt) then 58:       if (boxderivt) then
 59:          keep_angles = check_angles(box_params(4:6)) 59:          keep_angles = check_angles(box_params(4:6))
 60:          if (.not.keep_angles) then 60:          if (.not.keep_angles) then
 61:             ! reject the structure if unrealistic cell angles 61:             ! reject the structure if unrealistic cell angles
 62:             call reject(energy, g) 62:             call reject(energy, g)
 63:             return 63:             return
 64:          endif 64:          endif
 65:       endif 65:       endif
 66:  66: 
 67:       call build_H(H, H_grad, gtest) ! H matrix 67:       call build_H(H, H_grad, gtest) ! H matrix
430:                               ! total gradient wrt CoM coords for rigid body J1430:                               ! total gradient wrt CoM coords for rigid body J1
431:                               G(J3-2:J3) = G(J3-2:J3) + DVDR*RSS(:) + FRIJ(:)431:                               G(J3-2:J3) = G(J3-2:J3) + DVDR*RSS(:) + FRIJ(:)
432:                               ! total gradient wrt CoM coords for rigid body J2432:                               ! total gradient wrt CoM coords for rigid body J2
433:                               G(J4-2:J4) = G(J4-2:J4) - DVDR*RSS(:) - FRIJ(:)433:                               G(J4-2:J4) = G(J4-2:J4) - DVDR*RSS(:) - FRIJ(:)
434: 434: 
435:                               ! total gradient wrt AA coords for rigid body J1435:                               ! total gradient wrt AA coords for rigid body J1
436:                               G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPI(:) + TIJ(:)436:                               G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPI(:) + TIJ(:)
437:                               ! total gradient wrt AA coords for rigid body J2437:                               ! total gradient wrt AA coords for rigid body J2
438:                               G(J6-2:J6) = G(J6-2:J6) + DVDR*DRIJDPJ(:) + TJI(:)438:                               G(J6-2:J6) = G(J6-2:J6) + DVDR*DRIJDPJ(:) + TJI(:)
439: 439: 
440:                               ! compute gradients wrt lattice parameters440:                               ! gradients wrt cell parameters
441:                               if (boxderivt) then441:                               if (boxderivt) then
442: 442:                                  do idx = 1, 6
443:                               do idx = 1, 6443:                                     box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product((dvdr*rss(1:3) + frij(1:3)), matmul(H_grad(:,:,idx), rcomfrac(:)))
444:                                  box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product((dvdr*rss(1:3) + frij(1:3)), matmul(H_grad(:,:,idx), rcomfrac(:)))444:                                  enddo
445:                               enddo 
446:                               endif ! box derivatives445:                               endif ! box derivatives
447: 446: 
448:                            ENDIF ! gtest447:                            ENDIF ! gtest
 448:       
449:                         ENDIF ! within cutoff449:                         ENDIF ! within cutoff
450: 450: 
451:                      enddo ! n451:                      enddo ! n
452:                   enddo ! m452:                   enddo ! m
453:                enddo ! l453:                enddo ! l
454: 454: 
455:                ENDDO ! sites j455:                ENDDO ! sites j
456: 456: 
457:             ENDDO ! rigid bodies J457:             ENDDO ! rigid bodies J
458:  458:  
471:          ! loop over sites i471:          ! loop over sites i
472:          do i = 1, nsiteperbody(j1)472:          do i = 1, nsiteperbody(j1)
473:             j7 = maxsite*(j1-1) + i473:             j7 = maxsite*(j1-1) + i
474:             ei(:) = e(j7,:)474:             ei(:) = e(j7,:)
475: 475: 
476:             ! loop over sites j476:             ! loop over sites j
477:             do j = 1, nsiteperbody(j1)477:             do j = 1, nsiteperbody(j1)
478:                j8 = maxsite*(j1-1) + j478:                j8 = maxsite*(j1-1) + j
479:                ej(:) = e(j8,:)479:                ej(:) = e(j8,:)
480: 480: 
481:                ! get absolute displacement481:                ! site-site separation vector
482:                rr(:) = r(j7,:) - r(j8,:)482:                rr(:) = r(j7,:) - r(j8,:)
483:                ! convert to fractional483:                ! convert to fractional
484:                rrfrac(:) = matmul(H_inverse, rr(:))484:                rrfrac(:) = matmul(H_inverse, rr(:))
485: 485: 
486:                ! sum over lattice vectors486:                ! sum over lattice vectors
487:                do l = -newaldreal(1), newaldreal(1)487:                do l = -newaldreal(1), newaldreal(1)
488:                   do m = -newaldreal(2), newaldreal(2)488:                   do m = -newaldreal(2), newaldreal(2)
489:                      do n = -newaldreal(3), newaldreal(3)489:                      do n = -newaldreal(3), newaldreal(3)
490: 490: 
491:                      ! make sure not on same molecule491:                      ! if not on same rigid body
492:                      if (.not.(l.eq.0.and.m.eq.0.and.n.eq.0)) then492:                      if (.not.(l.eq.0.and.m.eq.0.and.n.eq.0)) then
493: 493: 
494:                         rssfrac(1) = rrfrac(1) + l494:                         rssfrac(1) = rrfrac(1) + l
495:                         rssfrac(2) = rrfrac(2) + m495:                         rssfrac(2) = rrfrac(2) + m
496:                         rssfrac(3) = rrfrac(3) + n496:                         rssfrac(3) = rrfrac(3) + n
497:                         ! convert back to absolute497:                         ! convert back to absolute
498:                         rss(:) = matmul(H, rssfrac(:))498:                         rss(:) = matmul(H, rssfrac(:))
499: 499: 
500:                         ! get COM displacement500:                         ! get vector between COM
501:                         if (gtest.and.boxderivt) then501:                         if (gtest.and.boxderivt) then
502:                            rcomfrac(1) = l502:                            rcomfrac(1) = l
503:                            rcomfrac(2) = m503:                            rcomfrac(2) = m
504:                            rcomfrac(3) = n504:                            rcomfrac(3) = n
505:                         endif505:                         endif
506: 506: 
507:                         r2 = dot_product(rss(:), rss(:))507:                         r2 = dot_product(rss(:), rss(:))
508:                         ! check within cutoff 
509:                         if (r2 < ewaldrealc2) then508:                         if (r2 < ewaldrealc2) then
510: 509: 
 510:                         ! absolute site-site distance
511:                         absrij = dsqrt(r2)511:                         absrij = dsqrt(r2)
512:                         nr(:) = rss(:)/absrij512:                         nr(:) = rss(:)/absrij
513:                         r2 = 1.d0/r2513:                         r2 = 1.d0/r2
514:                         r6 = r2*r2*r2514:                         r6 = r2*r2*r2
515: 515: 
516:                         ! CALCULATE DISPERSION DAMPING FACTOR516:                         ! CALCULATE DISPERSION DAMPING FACTOR
517: 517: 
518:                         ! initialize sum for the damping function and vertical shift518:                         ! initialize sum for the damping function and vertical shift
519:                         DMPFCT = 1.D0519:                         DMPFCT = 1.D0
520:                         DMPFCT_SHIFT = 1.D0520:                         DMPFCT_SHIFT = 1.D0
689:                         ENDIF689:                         ENDIF
690: 690: 
691: 691: 
692:                         IF (GTEST) THEN692:                         IF (GTEST) THEN
693: 693: 
694:                            ! total gradient wrt AA coords for rigid body J1694:                            ! total gradient wrt AA coords for rigid body J1
695:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPI(:) + TIJ(:)695:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPI(:) + TIJ(:)
696:                            ! total gradient wrt AA coords for rigid body J2696:                            ! total gradient wrt AA coords for rigid body J2
697:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPJ(:) + TJI(:)697:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPJ(:) + TJI(:)
698: 698: 
699:                            ! gradietn wrt lattice parameters699:                            ! gradients wrt cell parameters
700:                            if (boxderivt) then700:                            if (boxderivt) then
701:                            do idx = 1, 6701:                               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(:)))702:                                     box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product((dvdr*rss(1:3) + frij(1:3)), matmul(H_grad(:,:,idx), rcomfrac(:)))
703:                            enddo703:                               enddo
704:                            endif ! box derivatives704:                            endif ! box derivatives
705: 705: 
706:                         ENDIF ! gtest706:                         ENDIF ! gtest
707:                         endif ! central box707:                         endif ! central box
708:                     endif ! within cutoff708:                     endif ! within cutoff
709:                   enddo ! n709:                   enddo ! n
710:                enddo ! m710:                enddo ! m
711:             enddo ! l711:             enddo ! l
712:             enddo ! sites j712:             enddo ! sites j
713:          enddo ! sites i713:          enddo ! sites i
743:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XR, X)743:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XR, X)
744:          X(:) = XR(:)744:          X(:) = XR(:)
745:       ENDIF745:       ENDIF
746: 746: 
747:       ! add WCA-style repulsion to keep cell volume away from zero747:       ! 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)748:       if (boxderivt) call constrain_volume(v_fact, dv_fact, energy1, box_paramsgrad(4:6), gtest)
749: 749: 
750:       ! sum energies / gradients and convert to kJ/mol750:       ! sum energies / gradients and convert to kJ/mol
751:       ENERGY = (ENERGY1 + ENERGY2 + ENERGY3)*2625.499D0751:       ENERGY = (ENERGY1 + ENERGY2 + ENERGY3)*2625.499D0
752:       IF (GTEST) G(:) = (G(:) + G3(:))*2625.499D0752:       IF (GTEST) G(:) = (G(:) + G3(:))*2625.499D0
753:       if (gtest.and.boxderivt) box_paramsgrad(1:6) = box_paramsgrad(1:6)*2625.499D0753:       if (gtest) box_paramsgrad(1:6) = box_paramsgrad(1:6)*2625.499D0
754: 754: 
755:       END SUBROUTINE BENZGENRIGIDEWALD755:       END SUBROUTINE BENZGENRIGIDEWALD
756: 756: 
757: !     ----------------------------------------------------------------------------------------------757: !     ----------------------------------------------------------------------------------------------
758: !758: !
759: !      SUBROUTINE DEFPAHARIGID()759: !      SUBROUTINE DEFPAHARIGID()
760: !760: !
761: !      USE COMMONS, ONLY: RHOCC0, RHOCC10, RHOCC20,  RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, RHOCH10, RHOC20H, RHOCH20, &761: !      USE COMMONS, ONLY: RHOCC0, RHOCC10, RHOCC20,  RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, RHOCH10, RHOC20H, RHOCH20, &
762: !                         ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, CCKJ762: !                         ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, CCKJ
763: !763: !


r33133/benzgenrigid_ortho.f90 2017-08-07 15:30:36.209332643 +0100 r33132/benzgenrigid_ortho.f90 2017-08-07 15:30:38.413361964 +0100
 10:      &                   RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, RHOCH10, RHOC20H, & 10:      &                   RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, RHOCH10, RHOC20H, &
 11:      &                   RHOCH20, ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, & 11:      &                   RHOCH20, ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, &
 12:      &                   EWALDREALC, BOX_PARAMS, BOX_PARAMSGRAD 12:      &                   EWALDREALC, BOX_PARAMS, BOX_PARAMSGRAD
 13:  13: 
 14:       ! adapted to the genrigid framework 14:       ! adapted to the genrigid framework
 15:       USE GENRIGID, ONLY: NRIGIDBODY, ATOMRIGIDCOORDT, TRANSFORMCTORIGID, NSITEPERBODY, & 15:       USE GENRIGID, ONLY: NRIGIDBODY, ATOMRIGIDCOORDT, TRANSFORMCTORIGID, NSITEPERBODY, &
 16:      &                    MAXSITE, SITESRIGIDBODY, TRANSFORMRIGIDTOC, TRANSFORMGRAD 16:      &                    MAXSITE, SITESRIGIDBODY, TRANSFORMRIGIDTOC, TRANSFORMGRAD
 17:  17: 
 18:       ! use Ewald summation to compute electrostatics 18:       ! use Ewald summation to compute electrostatics
 19:       USE EWALD 19:       USE EWALD
 20:       ! box derivatives and fractional coordinates 
 21:       USE CARTDIST 20:       USE CARTDIST
 22:       USE BOX_DERIVATIVES 21:       USE BOX_DERIVATIVES
 23:  22: 
 24:       IMPLICIT NONE 23:       IMPLICIT NONE
 25:  24: 
 26:       INTEGER          :: I, J, K, J1, J2, J3, J4, J5, J6, J7, J8, OFFSET, FCT(6), L, M, N 25:       INTEGER          :: I, J, K, J1, J2, J3, J4, J5, J6, J7, J8, OFFSET, FCT(6), L, M, N
 27:       INTEGER          :: NEWALDREAL(3) 26:       INTEGER          :: NEWALDREAL(3)
 28:       DOUBLE PRECISION :: X(3*NATOMS) 27:       DOUBLE PRECISION :: X(3*NATOMS)
 29:       DOUBLE PRECISION, INTENT(OUT) :: G(3*NATOMS) 28:       DOUBLE PRECISION, INTENT(OUT) :: G(3*NATOMS)
 30:       DOUBLE PRECISION :: XR(3*NATOMS), XC(3*NATOMS), G3C(3*NATOMS), G3(3*NATOMS) 29:       DOUBLE PRECISION :: XR(3*NATOMS), XC(3*NATOMS), G3C(3*NATOMS), G3(3*NATOMS)
168:                      rrcom(:) = x(j3-2:j3) - x(j4-2:j4)167:                      rrcom(:) = x(j3-2:j3) - x(j4-2:j4)
169:                      ! minimum image convention168:                      ! minimum image convention
170:                      rrcommin(1) = rrcom(1) - box_params(1)*anint(rr(1)/box_params(1))169:                      rrcommin(1) = rrcom(1) - box_params(1)*anint(rr(1)/box_params(1))
171:                      rrcommin(2) = rrcom(2) - box_params(2)*anint(rr(2)/box_params(2))170:                      rrcommin(2) = rrcom(2) - box_params(2)*anint(rr(2)/box_params(2))
172:                      rrcommin(3) = rrcom(3) - box_params(3)*anint(rr(3)/box_params(3))171:                      rrcommin(3) = rrcom(3) - box_params(3)*anint(rr(3)/box_params(3))
173:                   endif172:                   endif
174: 173: 
175:                   ! sum over lattice vectors174:                   ! sum over lattice vectors
176:                   do l = -newaldreal(1), newaldreal(1)175:                   do l = -newaldreal(1), newaldreal(1)
177:                   rss(1) = rssmin(1) + box_params(1)*l176:                   rss(1) = rssmin(1) + box_params(1)*l
 177: 
178:                      do m = -newaldreal(2), newaldreal(2)178:                      do m = -newaldreal(2), newaldreal(2)
179:                      rss(2) = rssmin(2) + box_params(2)*m179:                      rss(2) = rssmin(2) + box_params(2)*m
180: 180: 
181:                         do n = -newaldreal(3), newaldreal(3)181:                         do n = -newaldreal(3), newaldreal(3)
182:                         rss(3) = rssmin(3) + box_params(3)*n182:                         rss(3) = rssmin(3) + box_params(3)*n
183: 183: 
184:                         ! get COM vector184:                         ! get COM vector
185:                         if (gtest.and.boxderivt) then185:                         if (gtest.and.boxderivt) then
186:                            rcom(1) = rrcommin(1) + box_params(1)*l186:                            rcom(1) = rrcommin(1) + box_params(1)*l
187:                            rcom(2) = rrcommin(2) + box_params(2)*m187:                            rcom(2) = rrcommin(2) + box_params(2)*m
439:                         rss(3) = rr(3) + box_params(3)*n439:                         rss(3) = rr(3) + box_params(3)*n
440: 440: 
441:                         ! get COM vector441:                         ! get COM vector
442:                         if (gtest.and.boxderivt) then442:                         if (gtest.and.boxderivt) then
443:                            rcom(1) = box_params(1)*l443:                            rcom(1) = box_params(1)*l
444:                            rcom(2) = box_params(2)*m444:                            rcom(2) = box_params(2)*m
445:                            rcom(3) = box_params(3)*n445:                            rcom(3) = box_params(3)*n
446:                         endif446:                         endif
447: 447: 
448:                         r2 = dot_product(rss(:), rss(:))448:                         r2 = dot_product(rss(:), rss(:))
449:                         ! check within cutoff 
450:                         if (r2 < ewaldrealc2) then449:                         if (r2 < ewaldrealc2) then
451: 450: 
 451:                         ! absolute site-site distance
452:                         absrij = dsqrt(r2)452:                         absrij = dsqrt(r2)
453:                         nr(:) = rss(:)/absrij453:                         nr(:) = rss(:)/absrij
454:                         r2 = 1.d0/r2454:                         r2 = 1.d0/r2
455:                         r6 = r2*r2*r2455:                         r6 = r2*r2*r2
456: 456: 
457:                         ! CALCULATE DISPERSION DAMPING FACTOR457:                         ! CALCULATE DISPERSION DAMPING FACTOR
458: 458: 
459:                         ! initialize sum for the damping function and vertical shift459:                         ! initialize sum for the damping function and vertical shift
460:                         DMPFCT = 1.D0460:                         DMPFCT = 1.D0
461:                         DMPFCT_SHIFT = 1.D0461:                         DMPFCT_SHIFT = 1.D0
622:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDR(:))622:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDR(:))
623:                               TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPI(:) &623:                               TIJ(:)  = ALPHACH*EXPFCT*(-DRIJDPI(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPI(:) &
624:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPI(:))624:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPI(:))
625:                               TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPJ(:) &625:                               TJI(:)  = ALPHACH*EXPFCT*(-DRIJDPJ(:)/ABSRIJ + (RHOCH10 + 3.D0*RHOCH20*COSTA)*DCADPJ(:) &
626:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPJ(:))626:                                       + (RHOC10H + 3.D0*RHOC20H*COSTB)*DCBDPJ(:))
627: 627: 
628:                            ENDIF628:                            ENDIF
629: 629: 
630:                         ENDIF630:                         ENDIF
631: 631: 
 632: 
632:                         IF (GTEST) THEN633:                         IF (GTEST) THEN
633: 634: 
634:                            ! total gradient wrt AA coords for rigid body J1635:                            ! total gradient wrt AA coords for rigid body J1
635:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPI(:) + TIJ(:)636:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPI(:) + TIJ(:)
636:                            ! total gradient wrt AA coords for rigid body J2637:                            ! total gradient wrt AA coords for rigid body J2
637:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPJ(:) + TJI(:)638:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPJ(:) + TJI(:)
638: 639: 
639:                            ! gradient wrt cell lengths640:                            ! gradient wrt cell lengths
640:                            if (boxderivt) box_paramsgrad(1:3) = box_paramsgrad(1:3) + (dvdr*rss(1:3)+frij(1:3))*rcom(1:3)/box_params(1:3)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)
641: 642: 


r33133/box_derivatives.f90 2017-08-07 15:30:36.429335570 +0100 r33132/box_derivatives.f90 2017-08-07 15:30:38.633364890 +0100
349: 349: 
350:           if (gtest) then350:           if (gtest) then
351:              ! add gradient contribution351:              ! 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(:)352:              grad_angles(:) = grad_angles(:) + 24.0d0*v_eps/v_sigma*((v_sigma/v)**7 - 2.0d0*(v_sigma/v)**13)*v_deriv(:)
353:           endif353:           endif
354:        endif354:        endif
355: 355: 
356:        end subroutine constrain_volume356:        end subroutine constrain_volume
357: 357: 
358: ! ---------------------------------------------------------------------------------358: ! ---------------------------------------------------------------------------------
 359: ! 
359: ! Rotates all rigid bodies after a step is taken in the cell parameters.360: ! Rotates all rigid bodies after a step is taken in the cell parameters.
360: ! VARIABLES361: ! VARIABLES
361: ! box_paramsold: unit cell lengths and angles from before the step was taken362: ! box_paramsold: unit cell lengths and angles from before the step was taken
362: ! TODO: this definitely is not working properly...363: ! TODO: not sure this is working properly (or that the equations are even valid...)
 364: 
 365: 
363: !       subroutine rotate_bodies(box_paramsold, xrfrac)366: !       subroutine rotate_bodies(box_paramsold, xrfrac)
364: !367: !
365: !       use genrigid, only: transformctorigid, sitesrigidbody, maxsite, nsiteperbody, inversematrix368: !       use genrigid, only: transformctorigid, sitesrigidbody, maxsite, nsiteperbody, inversematrix
366: !369: !
367: !       integer                         :: j1, j3, j5, j2, j4, j6370: !       integer                         :: j1, j3, j5, j2, j4, j6
368: !       double precision, intent(in)    :: box_paramsold(6)371: !       double precision, intent(in)    :: box_paramsold(6)
369: !       double precision, intent(inout) :: xrfrac(3*natoms)372: !       double precision, intent(inout) :: xrfrac(3*natoms)
370: !       double precision                :: H_old(3,3), H_grad(3,3,6), H_oldinverse(3,3), H_new(3,3)373: !       double precision                :: H_old(3,3), H_grad(3,3,6), H_oldinverse(3,3), H_new(3,3)
371: !       double precision                :: ri(3), p(3), xr(3*natoms), x(3*natoms), rot(3,3)374: !       double precision                :: ri(3), p(3), xr(3*natoms), x(3*natoms), rot(3,3)
372: !       double precision                :: rmi(3,3), drmi1(3,3), drmi2(3,3), drmi3(3,3), vol_new, vol_old375: !       double precision                :: rmi(3,3), drmi1(3,3), drmi2(3,3), drmi3(3,3), vol_new, vol_old


r33133/checkd.f90 2017-08-07 15:30:36.645338444 +0100 r33132/checkd.f90 2017-08-07 15:30:38.853367818 +0100
 13:       LOGICAL          :: GTEST, STEST, COMPON 13:       LOGICAL          :: GTEST, STEST, COMPON
 14:       DOUBLE PRECISION, PARAMETER :: ERRLIM = 1.D-04, DELX = 1.0D-6 14:       DOUBLE PRECISION, PARAMETER :: ERRLIM = 1.D-04, DELX = 1.0D-6
 15:       COMMON /CO/ COMPON 15:       COMMON /CO/ COMPON
 16:  16: 
 17:       ! dj337: allow for rigid bodies 17:       ! dj337: allow for rigid bodies
 18:       if (rigidinit) then 18:       if (rigidinit) then
 19:          dof = degfreedoms 19:          dof = degfreedoms
 20:       else 20:       else
 21:          dof = 3*natoms 21:          dof = 3*natoms
 22:       endif 22:       endif
  23: 
 23:       print *, 'DELX: ', DELX 24:       print *, 'DELX: ', DELX
 24:  25: 
 25: ! jwrm2> Turning compression on, if required 26: ! jwrm2> Turning compression on, if required
 26:       IF (COMPRESST .OR. PERCOLATET) COMPON = .TRUE. 27:       IF (COMPRESST .OR. PERCOLATET) COMPON = .TRUE.
 27:  28: 
 28: ! jwrm2> Converting GTHOMSON coordinates to polars 29: ! jwrm2> Converting GTHOMSON coordinates to polars
 29:       IF (GTHOMSONT) THEN 30:       IF (GTHOMSONT) THEN
 30:         CALL GTHOMSONCTOANG(X(1:3*NATOMS), TMPCOORDS(1:3*NATOMS), NATOMS) 31:         CALL GTHOMSONCTOANG(X(1:3*NATOMS), TMPCOORDS(1:3*NATOMS), NATOMS)
 31:         X(1:3*NATOMS) = TMPCOORDS(1:3*NATOMS) 32:         X(1:3*NATOMS) = TMPCOORDS(1:3*NATOMS)
 32:       END IF 33:       END IF
 47:          WRITE(*, *) IVRNO 48:          WRITE(*, *) IVRNO
 48:  49: 
 49:          ! dj337: make sure coordinates rigid body 50:          ! dj337: make sure coordinates rigid body
 50:          if (rigidinit.and.atomrigidcoordt) then 51:          if (rigidinit.and.atomrigidcoordt) then
 51:             call transformctorigid(x, tmpcoords) 52:             call transformctorigid(x, tmpcoords)
 52:             x(1:degfreedoms) = tmpcoords(1:degfreedoms) 53:             x(1:degfreedoms) = tmpcoords(1:degfreedoms)
 53:             x(degfreedoms+1:3*natoms) = 0.0d0 54:             x(degfreedoms+1:3*natoms) = 0.0d0
 54:             atomrigidcoordt = .false. 55:             atomrigidcoordt = .false.
 55:          endif 56:          endif
 56:  57: 
  58:          !print *, 'x: ', x(ivrno)
  59: 
 57:          GTEST = .FALSE. 60:          GTEST = .FALSE.
 58:          X(IVRNO) = X(IVRNO) - DELX 61:          X(IVRNO) = X(IVRNO) - DELX
  62:          !print *, 'xplus1: ', x(ivrno)
 59:          CALL POTENTIAL(X, G, FM, GTEST, STEST) 63:          CALL POTENTIAL(X, G, FM, GTEST, STEST)
  64:          !print *, 'xplus2: ', x(ivrno)
 60:  65: 
 61:          if (rigidinit.and.atomrigidcoordt) then 66:          if (rigidinit.and.atomrigidcoordt) then
 62:             call transformctorigid(x, tmpcoords) 67:             call transformctorigid(x, tmpcoords)
 63:             x(1:degfreedoms) = tmpcoords(1:degfreedoms) 68:             x(1:degfreedoms) = tmpcoords(1:degfreedoms)
 64:             x(degfreedoms+1:3*natoms) = 0.0d0 69:             x(degfreedoms+1:3*natoms) = 0.0d0
 65:             atomrigidcoordt = .false. 70:             atomrigidcoordt = .false.
 66:          endif 71:          endif
 67:  72: 
 68:          X(IVRNO) = X(IVRNO) + 2.D0*DELX 73:          X(IVRNO) = X(IVRNO) + 2.D0*DELX
  74:          !print *, 'xminus1: ', x(ivrno)
 69:          CALL POTENTIAL(X, G, FP, GTEST, STEST) 75:          CALL POTENTIAL(X, G, FP, GTEST, STEST)
  76:          !print *, 'xminus2: ', x(ivrno)
 70:       77:      
 71:          if (rigidinit.and.atomrigidcoordt) then 78:          if (rigidinit.and.atomrigidcoordt) then
 72:             call transformctorigid(x, tmpcoords) 79:             call transformctorigid(x, tmpcoords)
 73:             x(1:degfreedoms) = tmpcoords(1:degfreedoms) 80:             x(1:degfreedoms) = tmpcoords(1:degfreedoms)
 74:             x(degfreedoms+1:3*natoms) = 0.0d0 81:             x(degfreedoms+1:3*natoms) = 0.0d0
 75:             atomrigidcoordt = .false.  82:             atomrigidcoordt = .false. 
 76:          endif 83:          endif
 77:   84:  
 78:          GTEST = .TRUE. 85:          GTEST = .TRUE.
 79:          X(IVRNO) = X(IVRNO) - DELX 86:          X(IVRNO) = X(IVRNO) - DELX


r33133/ewald.f90 2017-08-07 15:30:36.861341316 +0100 r33132/ewald.f90 2017-08-07 15:30:39.077370796 +0100
 18:  18: 
 19: ! All equations for energy and gradient of Coulomb summation follow from: 19: ! All equations for energy and gradient of Coulomb summation follow from:
 20: ! Karasawa, N. and Goddard III, W. A. J. Phys. Chem., 93, 7320-7327 (1989). 20: ! Karasawa, N. and Goddard III, W. A. J. Phys. Chem., 93, 7320-7327 (1989).
 21:   21:  
 22: ! All input / output are in absolute Cartesian coordinates. 22: ! All input / output are in absolute Cartesian coordinates.
 23:  23: 
 24: ! Assuming all units for length, charge, and energy are in atomic units. 24: ! Assuming all units for length, charge, and energy are in atomic units.
 25:  25: 
 26: ! Works for either orthorhombic or triclinic unit cells. Computes energy gradient wrt 26: ! Works for either orthorhombic or triclinic unit cells. Computes energy gradient wrt
 27: ! cell parameters when BOXDERIVT keyword is true. 27: ! cell parameters when BOXDERIVT keyword is true.
 28:  
 29: ! NOTE: ERFC is the built-in complementary error function. If it is giving you an error 
 30: ! during compilation, try updating to a newer version of your compiler! 
 31: ! ----------------------------------------------------------------------------------- 28: ! -----------------------------------------------------------------------------------
 32:       subroutine ewaldsum(n, x, g, etot, gtest) 29:       subroutine ewaldsum(n, x, g, etot, gtest)
 33:  30: 
 34:       use cartdist, only: get_reciplatvec, build_H 31:       use cartdist, only: get_reciplatvec, build_H
 35:  32: 
 36:       implicit none 33:       implicit none
 37:  34: 
 38:       integer, intent(in)           :: n 35:       integer, intent(in)           :: n
 39:       integer                       :: newaldreal(3), newaldrecip(3) 36:       integer                       :: newaldreal(3), newaldrecip(3)
 40:       logical, intent(in)           :: gtest 37:       logical, intent(in)           :: gtest
 55:          ! orthorhombic unit cell 52:          ! orthorhombic unit cell
 56:          if (ortho) then 53:          if (ortho) then
 57:             ! determine number of lattice vectors to sum over 54:             ! determine number of lattice vectors to sum over
 58:             newaldreal(:) = floor(ewaldrealc/box_params(1:3) + 0.5d0) 55:             newaldreal(:) = floor(ewaldrealc/box_params(1:3) + 0.5d0)
 59:             ! compute real-space contribution to energy 56:             ! compute real-space contribution to energy
 60:             call coulombreal_ortho(x, newaldreal, etot) 57:             call coulombreal_ortho(x, newaldreal, etot)
 61:  58: 
 62:             ! determine number of reciprocal lattice vectors to sum over 59:             ! determine number of reciprocal lattice vectors to sum over
 63:             newaldrecip(:) = floor(ewaldrecipc*box_params(1:3)/(2.0d0*pi)) 60:             newaldrecip(:) = floor(ewaldrecipc*box_params(1:3)/(2.0d0*pi))
 64:             ! compute reciprocal-space contribution to energy 61:             ! compute reciprocal-space contribution to energy
 65:             !call coulombrecip_ortho(x, newaldrecip, etot) 62:             call coulombrecip_ortho(x, newaldrecip, etot)
 66:  63: 
 67:             if (gtest) then 64:             if (gtest) then
 68:                ! compute real-space contribution to gradient 65:                ! compute real-space contribution to gradient
 69:                call coulombrealgrad_ortho(x, newaldreal, g) 66:                call coulombrealgrad_ortho(x, newaldreal, g)
 70:  67: 
 71:                ! compute reciprocal-space contribution to gradient 68:                ! compute reciprocal-space contribution to gradient
 72:                !call coulombrecipgrad_ortho(x, newaldrecip, g) 69:                call coulombrecipgrad_ortho(x, newaldrecip, g)
 73:             endif 70:             endif
 74:          ! triclinic unit cell 71:          ! triclinic unit cell
 75:          else 72:          else
 76:             ! get reciprocal lattice vectors 73:             ! get reciprocal lattice vectors
 77:             call get_reciplatvec(reciplatvec, reciplatvec_grad, .false.) 74:             call get_reciplatvec(reciplatvec, reciplatvec_grad, .false.)
 78:             ! determine number of lattice vectors to sum over 75:             ! determine number of lattice vectors to sum over
 79:             newaldreal(1) = floor(ewaldrealc*dsqrt(sum(reciplatvec(1,:)**2))/(2.0d0*pi) + 0.5d0) 76:             newaldreal(1) = floor(ewaldrealc*dsqrt(sum(reciplatvec(1,:)**2))/(2.0d0*pi) + 0.5d0)
 80:             newaldreal(2) = floor(ewaldrealc*dsqrt(sum(reciplatvec(2,:)**2))/(2.0d0*pi) + 0.5d0) 77:             newaldreal(2) = floor(ewaldrealc*dsqrt(sum(reciplatvec(2,:)**2))/(2.0d0*pi) + 0.5d0)
 81:             newaldreal(3) = floor(ewaldrealc*dsqrt(sum(reciplatvec(3,:)**2))/(2.0d0*pi) + 0.5d0) 78:             newaldreal(3) = floor(ewaldrealc*dsqrt(sum(reciplatvec(3,:)**2))/(2.0d0*pi) + 0.5d0)
 82:             ! compute real-space contribution to energy 79:             ! compute real-space contribution to energy
851:                         q1 = stchrg(j1)848:                         q1 = stchrg(j1)
852:                         g(j3-2) = g(j3-2) - q1*q1*mul*rmin(1)849:                         g(j3-2) = g(j3-2) - q1*q1*mul*rmin(1)
853:                         g(j3-1) = g(j3-1) - q1*q1*mul*rmin(2)850:                         g(j3-1) = g(j3-1) - q1*q1*mul*rmin(2)
854:                         g(j3)   = g(j3)   - q1*q1*mul*rmin(3)851:                         g(j3)   = g(j3)   - q1*q1*mul*rmin(3)
855: 852: 
856:                         ! compute contribution to box derivatives853:                         ! compute contribution to box derivatives
857:                         if (boxderivt) then854:                         if (boxderivt) then
858:                            if (rigidinit) then855:                            if (rigidinit) then
859:                               box_paramsgrad(1:3) = box_paramsgrad(1:3) - q1*q1*mul*rmin(1:3)*rcom(1:3)/box_params(1:3)856:                               box_paramsgrad(1:3) = box_paramsgrad(1:3) - q1*q1*mul*rmin(1:3)*rcom(1:3)/box_params(1:3)
860:                            else ! not rigid bodies857:                            else ! not rigid bodies
861:                               box_paramsgrad(1:3) = box_paramsgrad(1:3) - q1*q1*mul*rmin(1:3)*rmin(1:3)/box_params(1:3)858:                               box_paramsgrad(1:3) = box_paramsgrad(1:3) + q1*q1*mul*rmin(1:3)*rmin(1:3)/box_params(1:3)
862:                            endif859:                            endif
863:                         endif 860:                         endif 
864: 861: 
865:                      enddo ! atoms862:                      enddo ! atoms
866:                   endif ! within cutoff863:                   endif ! within cutoff
867:                endif ! not in central box864:                endif ! not in central box
868:             enddo ! n865:             enddo ! n
869:          enddo ! m866:          enddo ! m
870:       enddo ! l867:       enddo ! l
871: 868: 


r33133/keywords.f 2017-08-07 15:30:37.089344349 +0100 r33132/keywords.f 2017-08-07 15:30:39.309373882 +0100
3654:          IF (NITEMS.GT.4) CALL READF(RSPEED)3654:          IF (NITEMS.GT.4) CALL READF(RSPEED)
3655: 3655: 
3656:          ewaldalpha = 5.6d0/minval(box3d)3656:          ewaldalpha = 5.6d0/minval(box3d)
3657: 3657: 
3658:          NRBSITES = 123658:          NRBSITES = 12
3659:          ALLOCATE(SITE(NRBSITES,3))3659:          ALLOCATE(SITE(NRBSITES,3))
3660:          ALLOCATE(RBSTLA(NRBSITES,3))3660:          ALLOCATE(RBSTLA(NRBSITES,3))
3661: 3661: 
3662:          CALL DEFPAHARIGID()3662:          CALL DEFPAHARIGID()
3663:          NCARBON = 63663:          NCARBON = 6
 3664:          !nrigidbody = 32
 3665:          !CALL DEFBENZENERIGID()
3664: 3666: 
3665: ! ----------------------------------------------------------------------------------------3667: ! ----------------------------------------------------------------------------------------
3666: 3668: 
3667:       ELSE IF (WORD.EQ.'FAL') THEN3669:       ELSE IF (WORD.EQ.'FAL') THEN
3668:          FAL=.TRUE.3670:          FAL=.TRUE.
3669: 3671: 
3670:       ELSE IF (WORD == 'FEBH') THEN3672:       ELSE IF (WORD == 'FEBH') THEN
3671:          CALL READF(FETEMP)3673:          CALL READF(FETEMP)
3672:          FEBHT = .TRUE.3674:          FEBHT = .TRUE.
3673:          FE_FILE_UNIT = GETUNIT()3675:          FE_FILE_UNIT = GETUNIT()


r33133/main.F 2017-08-07 15:30:37.309347277 +0100 r33132/main.F 2017-08-07 15:30:39.529376809 +0100
200:               STOP                200:               STOP                
201:           END IF201:           END IF
202:       END IF202:       END IF
203: 203: 
204: ! dj337: allocate charges for benzenes204: ! dj337: allocate charges for benzenes
205:       if (benzrigidewaldt.or.ewaldt) then205:       if (benzrigidewaldt.or.ewaldt) then
206:          ALLOCATE(STCHRG(NRIGIDBODY*NRBSITES))206:          ALLOCATE(STCHRG(NRIGIDBODY*NRBSITES))
207:          CALL DEFBENZENERIGIDEWALD()207:          CALL DEFBENZENERIGIDEWALD()
208:       endif208:       endif
209: 209: 
 210: 
210:       IF (MULTIPOTT) THEN211:       IF (MULTIPOTT) THEN
211:           CALL MULTIPOT_INITIALISE()212:           CALL MULTIPOT_INITIALISE()
212:       ENDIF213:       ENDIF
213: 214: 
214:       IF (CUDAT) THEN215:       IF (CUDAT) THEN
215:          IF ((CUDAPOT .EQ. 'A') .AND. (.NOT. AMBER12T)) THEN216:          IF ((CUDAPOT .EQ. 'A') .AND. (.NOT. AMBER12T)) THEN
216:             WRITE(MYUNIT,'(A)') " main> The AMBER12 keyword must be used with 'CUDA A'. "217:             WRITE(MYUNIT,'(A)') " main> The AMBER12 keyword must be used with 'CUDA A'. "
217:             STOP218:             STOP
218:          END IF219:          END IF
219:          IF (DEBUG) THEN220:          IF (DEBUG) THEN
578: 579: 
579:       CALL FLUSH(6)580:       CALL FLUSH(6)
580: 581: 
581:       IF (PROJIT) THEN582:       IF (PROJIT) THEN
582:          ALLOCATE(TMPCOORDS(3*NATOMSALLOC))583:          ALLOCATE(TMPCOORDS(3*NATOMSALLOC))
583:          TMPCOORDS(1:3*NATOMSALLOC)=COORDS(1:3*NATOMSALLOC,1)584:          TMPCOORDS(1:3*NATOMSALLOC)=COORDS(1:3*NATOMSALLOC,1)
584:          CALL PROJIINIT(TMPCOORDS,NATOMSALLOC)585:          CALL PROJIINIT(TMPCOORDS,NATOMSALLOC)
585:          DEALLOCATE(TMPCOORDS)586:          DEALLOCATE(TMPCOORDS)
586:       ENDIF587:       ENDIF
587: 588: 
 589: 
588:       IF(VGW) THEN590:       IF(VGW) THEN
589:         CALL INITIALIZE_VGWSP(NATOMSALLOC,LJSIGMA,LJEPSILON,TAUMAX,CPFACTORSG,CPS,VGWTOL)591:         CALL INITIALIZE_VGWSP(NATOMSALLOC,LJSIGMA,LJEPSILON,TAUMAX,CPFACTORSG,CPS,VGWTOL)
590:         CALL INITIALIZE_VGW(NATOMSALLOC,LJSIGMA,LJEPSILON,TAUMAXFULL,CPFACTORFG,CPF)592:         CALL INITIALIZE_VGW(NATOMSALLOC,LJSIGMA,LJEPSILON,TAUMAXFULL,CPFACTORFG,CPF)
591:       ENDIF593:       ENDIF
592: 594: 
593:       IF (SEEDT) THEN595:       IF (SEEDT) THEN
594:          CALL GSEED596:          CALL GSEED
595:       ELSE597:       ELSE
596:          IF ((.NOT.FIELDT).AND.CENT) THEN598:          IF ((.NOT.FIELDT).AND.CENT) THEN
597:             DO J1=1,NPAR599:             DO J1=1,NPAR


r33133/mc.F 2017-08-07 15:30:37.533350255 +0100 r33132/mc.F 2017-08-07 15:30:39.753379789 +0100
2080:                            HBONDGROUPS(HBONDMATCHN,:,:)=HBONDMAT2080:                            HBONDGROUPS(HBONDMATCHN,:,:)=HBONDMAT
2081: ! csw34> Dump an updates group*.mat file for the group2081: ! csw34> Dump an updates group*.mat file for the group
2082:                            WRITE(NHBONDGROUPSCHAR,*) HBONDMATCHN2082:                            WRITE(NHBONDGROUPSCHAR,*) HBONDMATCHN
2083:                            HBONDGROUPNAME='group'//TRIM(ADJUSTL(NHBONDGROUPSCHAR))2083:                            HBONDGROUPNAME='group'//TRIM(ADJUSTL(NHBONDGROUPSCHAR))
2084:                            CALL WRITEHBONDMATRIX(HBONDMAT(:,:),HBONDNRES,HBONDGROUPNAME)2084:                            CALL WRITEHBONDMATRIX(HBONDMAT(:,:),HBONDNRES,HBONDGROUPNAME)
2085:                         ENDIF2085:                         ENDIF
2086:                      ENDIF2086:                      ENDIF
2087:                   ENDIF2087:                   ENDIF
2088: ! END OF HBONDMATRIX2088: ! END OF HBONDMATRIX
2089:                   COORDSO(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)2089:                   COORDSO(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)
2090:                   ! dj337: update saved box parameters2090:                   if (boxderivt) box_paramso(1:6) = box_params(1:6) ! dj337: update saved box parameters
2091:                   if (boxderivt) box_paramso(1:6) = box_params(1:6) 
2092:                   LABELSO(1:NATOMS,JP)=LABELS(1:NATOMS,JP) ! <ds6562091:                   LABELSO(1:NATOMS,JP)=LABELS(1:NATOMS,JP) ! <ds656
2093:                   VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)2092:                   VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)
2094:                ELSE2093:                ELSE
2095:                   NFAIL(JP)=NFAIL(JP)+12094:                   NFAIL(JP)=NFAIL(JP)+1
2096:                   REJSTREAK=REJSTREAK+12095:                   REJSTREAK=REJSTREAK+1
2097:                   MARKOVSTEP=.FALSE.2096:                   MARKOVSTEP=.FALSE.
2098:                   CALL MYRESET(JP,NATOMS,NPAR,NSEED)2097:                   CALL MYRESET(JP,NATOMS,NPAR,NSEED)
2099:                   IF (DEBUG) THEN2098:                   IF (DEBUG) THEN
2100:                      WRITE(MYUNIT,36) JP,RANDOM,POTEL,EPREV(JP),NSUCCESS(JP),NFAIL(JP)2099:                      WRITE(MYUNIT,36) JP,RANDOM,POTEL,EPREV(JP),NSUCCESS(JP),NFAIL(JP)
2101: 36                   FORMAT('JP,RAN,POTEL,EPREV,NSUC,NFAIL=',I2,3F15.7,2I6,' REJ')2100: 36                   FORMAT('JP,RAN,POTEL,EPREV,NSUC,NFAIL=',I2,3F15.7,2I6,' REJ')


r33133/potential.f90 2017-08-07 15:30:37.753353186 +0100 r33132/potential.f90 2017-08-07 15:30:39.977382766 +0100
343: !         EMINUS=EDUM1+EDUM2343: !         EMINUS=EDUM1+EDUM2
344: !344: !
345: !         X(J1)=X(J1)+DIFF345: !         X(J1)=X(J1)+DIFF
346: !         GRAD(J1)=(EPLUS-EMINUS)/(2.0D0*DIFF)346: !         GRAD(J1)=(EPLUS-EMINUS)/(2.0D0*DIFF)
347: !         WRITE(*,'(I5, 4F20.10)') J1, GRAD(J1), EPLUS, EMINUS, ABS(EPLUS-EMINUS)347: !         WRITE(*,'(I5, 4F20.10)') J1, GRAD(J1), EPLUS, EMINUS, ABS(EPLUS-EMINUS)
348: !     END DO348: !     END DO
349: !349: !
350: 350: 
351: ! dj337: pahagenrigid for benzene with Ewald351: ! dj337: pahagenrigid for benzene with Ewald
352:    ELSE IF (BENZRIGIDEWALDT) THEN352:    ELSE IF (BENZRIGIDEWALDT) THEN
353:       IF (ORTHO) THEN353:       if (ortho) then
354:          CALL BENZGENRIGIDEWALD_ORTHO(X, GRAD, EREAL, GRADT)354:          call benzgenrigidewald_ortho(x, grad, ereal, gradt)
355:       ELSE355:       else
356:          CALL BENZGENRIGIDEWALD(X, GRAD, EREAL, GRADT)356:          CALL BENZGENRIGIDEWALD(X, GRAD, EREAL, GRADT)
357:       ENDIF357:       endif
358:       if (rigidinit .and. (.not.atomrigidcoordt)) then358:       if (rigidinit .and. (.not.atomrigidcoordt)) then
359:          xrigidcoords(1:degfreedoms) = x(1:degfreedoms)359:          xrigidcoords(1:degfreedoms) = x(1:degfreedoms)
360:          xrigidgrad(1:degfreedoms) = grad(1:degfreedoms)360:          xrigidgrad(1:degfreedoms) = grad(1:degfreedoms)
361:          call transformrigidtoc(1, nrigidbody, xcoords, xrigidcoords)361:          call transformrigidtoc(1, nrigidbody, xcoords, xrigidcoords)
362:       endif362:       endif
363: ! check analytical and numerical gradients363: ! check analytical and numerical gradients
364: !      DIFF=1.0D-6364: !      DIFF=1.0D-6
365: !      WRITE(*, *) 'analytic and numerical gradients:'365: !      WRITE(*, *) 'analytic and numerical gradients:'
366: !      print *, 'diff: ', diff366: !      print *, 'diff: ', diff
367: !      IF (ATOMRIGIDCOORDT) THEN ! if input is cartesian367: !      IF (ATOMRIGIDCOORDT) THEN ! if input is cartesian


r33133/saveit.f 2017-08-07 15:30:37.973356108 +0100 r33132/saveit.f 2017-08-07 15:30:40.197385698 +0100
 49: C 49: C
 50: C  These are probably the same - but just to make sure we save the  50: C  These are probably the same - but just to make sure we save the 
 51: C  lowest. 51: C  lowest.
 52: C 52: C
 53:             IF (EREAL.LT.QMIN(J1)) THEN 53:             IF (EREAL.LT.QMIN(J1)) THEN
 54:                QMINNATOMS(J1)=NATOMS 54:                QMINNATOMS(J1)=NATOMS
 55:                QMIN(J1)=EREAL 55:                QMIN(J1)=EREAL
 56:                DO J2=1,3*NATOMS 56:                DO J2=1,3*NATOMS
 57:                   QMINP(J1,J2)=P(J2) 57:                   QMINP(J1,J2)=P(J2)
 58:                ENDDO 58:                ENDDO
 59:                ! dj337: save the box parameters 59:                ! dj337: save box parameters
 60:                if (boxderivt) boxq(j1,1:6) = box_params(1:6) 60:                if (boxderivt) boxq(j1,1:6) = box_params(1:6)
 61:             ENDIF 61:             ENDIF
 62:             GOTO 10 62:             GOTO 10
 63:          ENDIF 63:          ENDIF
 64:          IF (EREAL.LT.QMIN(J1)) THEN 64:          IF (EREAL.LT.QMIN(J1)) THEN
 65:  65: 
 66:             J2=NSAVE 66:             J2=NSAVE
 67: 20          CONTINUE 67: 20          CONTINUE
 68:             IF (NSAVE.GT.1) THEN 68:             IF (NSAVE.GT.1) THEN
 69:                QMIN(J2)=QMIN(J2-1) 69:                QMIN(J2)=QMIN(J2-1)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0