hdiff output

r33131/benzgenrigid.f90 2017-08-07 15:30:26.657205572 +0100 r33130/benzgenrigid.f90 2017-08-07 15:30:28.697232711 +0100
 17:       USE EWALD 17:       USE EWALD
 18:       USE CARTDIST 18:       USE CARTDIST
 19:       USE BOX_DERIVATIVES 19:       USE BOX_DERIVATIVES
 20:  20: 
 21:       IMPLICIT NONE 21:       IMPLICIT NONE
 22:  22: 
 23:       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
 24:       INTEGER          :: NEWALDREAL(3) 24:       INTEGER          :: NEWALDREAL(3)
 25:       DOUBLE PRECISION :: X(3*NATOMS) 25:       DOUBLE PRECISION :: X(3*NATOMS)
 26:       DOUBLE PRECISION, INTENT(OUT) :: G(3*NATOMS) 26:       DOUBLE PRECISION, INTENT(OUT) :: G(3*NATOMS)
 27:       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), graddum(3*natoms)
 28:       DOUBLE PRECISION, INTENT(OUT) :: ENERGY 28:       DOUBLE PRECISION, INTENT(OUT) :: ENERGY
 29:       DOUBLE PRECISION :: R2, R6, ABSRIJ, DVDR, ENERGY1, ENERGY2, ENERGY3 29:       DOUBLE PRECISION :: R2, R6, ABSRIJ, DVDR, ENERGY1, ENERGY2, ENERGY3, diff, eplus, eminus
 30:       DOUBLE PRECISION :: DMPFCT_SHIFT, EXPFCT_SHIFT, VSHIFT1, VSHIFT2, EWALDREALC2 30:       DOUBLE PRECISION :: DMPFCT_SHIFT, EXPFCT_SHIFT, VSHIFT1, VSHIFT2, EWALDREALC2, RCOMMIN(3), RCOM(3)
 31:       DOUBLE PRECISION :: RI(3), RR(3), RSS(3), NR(3), P(3), EI(3), EJ(3), FRIJ(3), TIJ(3), TJI(3)  31:       DOUBLE PRECISION :: RI(3), RR(3), RSS(3), RSSMIN(3), NR(3), P(3), EI(3), EJ(3), FRIJ(3), TIJ(3), TJI(3) 
 32:       DOUBLE PRECISION :: R(MAXSITE*NRIGIDBODY,3), E(3*MAXSITE*NRIGIDBODY,3), xdum(3*natoms) 32:       DOUBLE PRECISION :: R(MAXSITE*NRIGIDBODY,3), E(3*MAXSITE*NRIGIDBODY,3), xdum(3*natoms)
 33:       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)
 34:       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)
 35:       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)
 36:       DOUBLE PRECISION :: RHOCC, RHOHH, RHOCH, COSTA, COSTB, DMPFCT, DDMPDR, EXPFCT  36:       DOUBLE PRECISION :: RHOCC, RHOHH, RHOCH, COSTA, COSTB, DMPFCT, DDMPDR, EXPFCT 
 37:       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)
 38:       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)
 39:       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)
 40:       double precision :: reciplatvec(3,3), reciplatvec_grad(3,3,6) 40:       double precision :: H_mat(3,3), rmatrix(3), reciplatvec(3,3), reciplatvec_grad(3,3,6)
 41:       DOUBLE PRECISION, PARAMETER :: B = 1.6485D0 41:       DOUBLE PRECISION, PARAMETER :: B = 1.6485D0
 42:       double precision, parameter :: pi = 3.141592654d0 42:       double precision, parameter :: pi = 3.141592654d0
 43:       integer, parameter          :: image_cutoff = 5 43:       integer, parameter          :: image_cutoff = 5
 44:       LOGICAL          :: GTEST, keep_angles 44:       LOGICAL          :: GTEST, keep_angles
 45:  45: 
  46:       !print *, 'box_params: ', box_params(1:6)
 46:  47: 
 47:       if (boxderivt) then 48:       keep_angles = check_angles(box_params(4:6))
 48:          keep_angles = check_angles(box_params(4:6)) 49: 
 49:          if (.not.keep_angles) then 50:       if (.not.keep_angles) then
 50:             print *, 'rejecting1' 51:          !print *, 'rejecting'
 51:             call reject(energy, g) 52:          call reject(energy, g)
 52:             return 53:          return
 53:          endif 
 54:       endif 54:       endif
 55:  55: 
 56:       call build_H(H, H_grad, gtest) 56:       call build_H(H, H_grad, gtest)
 57:       call inversematrix(H, H_inverse) 57:       call inversematrix(H, H_inverse)
 58:       call get_volume(vol) 58:       call get_volume(vol)
 59:       call get_reciplatvec(reciplatvec,reciplatvec_grad, .false.) 59:       call get_reciplatvec(reciplatvec,reciplatvec_grad, .false.)
 60:  60: 
 61:       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)
 62:       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)
 63:       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) 64:       !print *, 'newaldreal: ', newaldreal(:3)
 65:  65: 
 66:       if (boxderivt) then 66:       if (.not. all(newaldreal.le.image_cutoff)) then
 67:          if (.not. all(newaldreal.le.image_cutoff)) then 67:          !print *, 'rejecting'
 68:             print *, 'rejecting2' 68:          call reject(energy, g)
 69:             call reject(energy, g) 69:          return
 70:             return 
 71:          endif 
 72:       endif 70:       endif
 73:  71: 
 74:       ! factorials 72:       ! factorials
 75:       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
 76:       ! initialize energy values 74:       ! initialize energy values
 77:       ! energy1 is due to short-range anisotropic interactions 75:       ! energy1 is due to short-range anisotropic interactions
 78:       ! energy2 is due to damped dispersion 76:       ! energy2 is due to damped dispersion
 79:       ! energy3 is due to long-range electrostatics (computed using Ewald) 77:       ! energy3 is due to long-range electrostatics (computed using Ewald)
 80:       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
 81:  79: 
104: 102: 
105:       !print *, 'coords benzrigid   : '103:       !print *, 'coords benzrigid   : '
106:       !do j1 = 1, 2*nrigidbody104:       !do j1 = 1, 2*nrigidbody
107:       !   print *, x(3*j1-2:3*j1)105:       !   print *, x(3*j1-2:3*j1)
108:       !enddo106:       !enddo
109: 107: 
110:       !call build_H(H, H_grad, gtest)108:       !call build_H(H, H_grad, gtest)
111:       !call inversematrix(H, H_inverse)109:       !call inversematrix(H, H_inverse)
112:       !call get_volume(vol)110:       !call get_volume(vol)
113:       !call get_reciplatvec(reciplatvec,reciplatvec_grad, .false.)111:       !call get_reciplatvec(reciplatvec,reciplatvec_grad, .false.)
114:       if (boxderivt) then112:       ! compute v factor
115:          ! compute v factor113:       c(:) = dcos(box_params(4:6))
116:          c(:) = dcos(box_params(4:6))114:       s(:) = dsin(box_params(4:6))
117:          s(:) = dsin(box_params(4:6))115:       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))116:       dv_fact(1) = s(1)*(c(1) - c(2)*c(3))/v_fact
119:          dv_fact(1) = s(1)*(c(1) - c(2)*c(3))/v_fact117:       dv_fact(2) = s(2)*(c(2) - c(1)*c(3))/v_fact
120:          dv_fact(2) = s(2)*(c(2) - c(1)*c(3))/v_fact118:       dv_fact(3) = s(3)*(c(3) - c(1)*c(2))/v_fact
121:          dv_fact(3) = s(3)*(c(3) - c(1)*c(2))/v_fact 
122:       endif 
123: 119: 
124:       EWALDREALC2 = EWALDREALC**2120:       EWALDREALC2 = EWALDREALC**2
125: 121: 
126:       !print *, 'reciplatvec: ', reciplatvec(:3,:3)122:       !print *, 'reciplatvec: ', reciplatvec(:3,:3)
127: 123: 
128:       ! OFFSET is number of CoM coords (3*NRIGIDBODY)124:       ! OFFSET is number of CoM coords (3*NRIGIDBODY)
129:       OFFSET     = 3*NRIGIDBODY125:       OFFSET     = 3*NRIGIDBODY
130: 126: 
131:       ! Computing Cartesian coordinates for the system.  127:       ! Computing Cartesian coordinates for the system.  
132:       DO J1 = 1, NRIGIDBODY128:       DO J1 = 1, NRIGIDBODY
235:                         rssfrac(3) = rssfracmin(3) + n231:                         rssfrac(3) = rssfracmin(3) + n
236: 232: 
237:                         ! convert to absolute coordinates233:                         ! convert to absolute coordinates
238:                         rss(:) = matmul(H, rssfrac(:))234:                         rss(:) = matmul(H, rssfrac(:))
239: 235: 
240:                         ! get COM vector236:                         ! get COM vector
241:                         if (gtest.and.boxderivt) then237:                         if (gtest.and.boxderivt) then
242:                            rcomfrac(1) = rcomfracmin(1) + l238:                            rcomfrac(1) = rcomfracmin(1) + l
243:                            rcomfrac(2) = rcomfracmin(2) + m239:                            rcomfrac(2) = rcomfracmin(2) + m
244:                            rcomfrac(3) = rcomfracmin(3) + n240:                            rcomfrac(3) = rcomfracmin(3) + n
 241:                            ! convert to absolute coordinates
 242:                            rcom(:) = matmul(H, rcomfrac(:))
245:                         endif243:                         endif
246:                      244:                      
247:                         R2     = DOT_PRODUCT(RSS(:),RSS(:))245:                         R2     = DOT_PRODUCT(RSS(:),RSS(:))
248:                         !print *, 'r2: ', r2246:                         !print *, 'r2: ', r2
249:                         ! check if distance within cutoff247:                         ! check if distance within cutoff
250:                         IF (R2 < EWALDREALC2) THEN248:                         IF (R2 < EWALDREALC2) THEN
251:                            !print *, j7, j8249:                            !print *, j7, j8
252:                            !print *, 'r   : ', rss(:3)250:                            !print *, 'r   : ', rss(:3)
253:                            !print *, 'rmin: ', rssmin(1:3)251:                            !print *, 'rmin: ', rssmin(1:3)
254:                            !print *, 'rcom: ', rcom(1:3)252:                            !print *, 'rcom: ', rcom(1:3)
461:                               ! total gradient wrt AA coords for rigid body J2459:                               ! total gradient wrt AA coords for rigid body J2
462:                               G(J6-2:J6) = G(J6-2:J6) + DVDR*DRIJDPJ(:) + TJI(:)460:                               G(J6-2:J6) = G(J6-2:J6) + DVDR*DRIJDPJ(:) + TJI(:)
463:                               !g(j6-2:j6) = g(j6-2:j6) + tji(:)461:                               !g(j6-2:j6) = g(j6-2:j6) + tji(:)
464: 462: 
465:                               ! TODO: this is if orientation of rigid body depends on cell parameters463:                               ! 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)464:                               !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(:)))465:                               !rmatrix(:) = matmul(H_inverse, vol*(rss(:)-rcom(:)))
468:                               !box_paramsgrad(4) = box_paramsgrad(4) + dvdr*dot_product(rss(1:3), matmul(H_mat, rmatrix(:)))466:                               !box_paramsgrad(4) = box_paramsgrad(4) + dvdr*dot_product(rss(1:3), matmul(H_mat, rmatrix(:)))
469: 467: 
470:                               !box_paramsgrad(1) = box_paramsgrad(1) + dvdr*dot_product(rss(1:3), matmul(H_mat/vol, rmatrix(:)))468:                               !box_paramsgrad(1) = box_paramsgrad(1) + dvdr*dot_product(rss(1:3), matmul(H_mat/vol, rmatrix(:)))
471:                               if (boxderivt) then469: 
472:                                  do idx = 1, 6470:                               do idx = 1, 6
473:                                     box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product((dvdr*rss(1:3) + frij(1:3)), matmul(H_grad(:,:,idx), rcomfrac(:)))471:                                  box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product((dvdr*rss(1:3) + frij(1:3)), matmul(H_grad(:,:,idx), rcomfrac(:)))
474:                                  enddo472:                               enddo
475:                               endif 
476: 473: 
477:                            ENDIF474:                            ENDIF
478:       475:       
479:                         ENDIF476:                         ENDIF
480: 477: 
481:                      enddo478:                      enddo
482:                   enddo479:                   enddo
483:                enddo480:                enddo
484: 481: 
485:                ENDDO482:                ENDDO
528:                         rssfrac(2) = rrfrac(2) + m525:                         rssfrac(2) = rrfrac(2) + m
529:                         rssfrac(3) = rrfrac(3) + n526:                         rssfrac(3) = rrfrac(3) + n
530:                         !print *, l, m, n527:                         !print *, l, m, n
531:                         !print *, 'rssfrac: ', rssfrac(:3)528:                         !print *, 'rssfrac: ', rssfrac(:3)
532:                         rss(:) = matmul(H, rssfrac(:))529:                         rss(:) = matmul(H, rssfrac(:))
533: 530: 
534:                         if (gtest.and.boxderivt) then531:                         if (gtest.and.boxderivt) then
535:                            rcomfrac(1) = l532:                            rcomfrac(1) = l
536:                            rcomfrac(2) = m533:                            rcomfrac(2) = m
537:                            rcomfrac(3) = n534:                            rcomfrac(3) = n
 535:                            rcom(:) = matmul(H, rcomfrac(:))
538:                         endif536:                         endif
539: 537: 
540:                         r2 = dot_product(rss(:), rss(:))538:                         r2 = dot_product(rss(:), rss(:))
541:                         !print *, 'rssmin2: ', rss(1:3)539:                         !print *, 'rssmin2: ', rss(1:3)
542:                         if (r2 < ewaldrealc2) then540:                         if (r2 < ewaldrealc2) then
543: 541: 
544:                         !print *, 'rssmin3: ', rss(1:3)542:                         !print *, 'rssmin3: ', rss(1:3)
545:                         !print *, 'r2    : ', r2543:                         !print *, 'r2    : ', r2
546:                         !print *, 'rcom  : ', rcom(1:3)544:                         !print *, 'rcom  : ', rcom(1:3)
547:                         absrij = dsqrt(r2)545:                         absrij = dsqrt(r2)
557:                         DMPFCT = 1.D0555:                         DMPFCT = 1.D0
558:                         DMPFCT_SHIFT = 1.D0556:                         DMPFCT_SHIFT = 1.D0
559:                         ! initialize sum for the derivative of damping function557:                         ! initialize sum for the derivative of damping function
560:                         DDMPDR = B558:                         DDMPDR = B
561: 559: 
562:                         ! calculate sums560:                         ! calculate sums
563:                         DO K = 1, 6561:                         DO K = 1, 6
564: 562: 
565:                            DMPFCT = DMPFCT + (B*ABSRIJ)**K/FLOAT(FCT(K))563:                            DMPFCT = DMPFCT + (B*ABSRIJ)**K/FLOAT(FCT(K))
566:                            DMPFCT_SHIFT = DMPFCT_SHIFT + (B*EWALDREALC)**K/FLOAT(FCT(K))564:                            DMPFCT_SHIFT = DMPFCT_SHIFT + (B*EWALDREALC)**K/FLOAT(FCT(K))
567:                            IF (K > 1) DDMPDR = DDMPDR + (B**K)*(ABSRIJ)**(K-1)/FLOAT(FCT(K-1))565:                            !IF (K > 1) DDMPDR = DDMPDR + (B**K)*(ABSRIJ)**(K-1)/FLOAT(FCT(K-1))
568: 566: 
569:                         END DO567:                         END DO
570: 568: 
571:                         EXPFCT = DEXP(-B*ABSRIJ)569:                         EXPFCT = DEXP(-B*ABSRIJ)
572:                         EXPFCT_SHIFT = DEXP(-B*EWALDREALC)570:                         EXPFCT_SHIFT = DEXP(-B*EWALDREALC)
573:                         ! DDMPDR is derivative of damping function with factor 1/Rab571:                         ! DDMPDR is derivative of damping function with factor 1/Rab
574:                         DDMPDR = (B*EXPFCT*DMPFCT - EXPFCT*DDMPDR)/ABSRIJ572:                         DDMPDR = (B*EXPFCT*DMPFCT - EXPFCT*DDMPDR)/ABSRIJ
575:                         ! DMPFCT is damping function573:                         ! DMPFCT is damping function
576:                         DMPFCT = 1.D0 - EXPFCT*DMPFCT574:                         DMPFCT = 1.D0 - EXPFCT*DMPFCT
577:                         ! DMPFCT_SHIFT is vertical shift for damping function575:                         ! DMPFCT_SHIFT is vertical shift for damping function
748:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPJ(:) + TJI(:)746:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPJ(:) + TJI(:)
749:                            !g(j5-2:j5) = g(j5-2:j5) + tji(:)747:                            !g(j5-2:j5) = g(j5-2:j5) + tji(:)
750:                            !print *, 'dispersion:'748:                            !print *, 'dispersion:'
751:                            !print *, dvdr*drijdpi(:3)749:                            !print *, dvdr*drijdpi(:3)
752:                            !print *, dvdr*drijdpj(:3)750:                            !print *, dvdr*drijdpj(:3)
753:                            !print *, 'energy2: ', energy2751:                            !print *, 'energy2: ', energy2
754:                            !print *, 'anisotropic:'752:                            !print *, 'anisotropic:'
755:                            !print *, tij(:3)753:                            !print *, tij(:3)
756:                            !print *, tji(:3)754:                            !print *, tji(:3)
757: 755: 
758:                            if (boxderivt) then756:                            do idx = 1, 6
759:                               do idx = 1, 6757:                                  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(:)))758:                                  !box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product(dvdr*rss(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(:)))759:                            enddo
762:                               enddo 
763:                            endif 
764: 760: 
765:                         ENDIF761:                         ENDIF
766:                         endif762:                         endif
767:                     endif763:                     endif
768:                   enddo764:                   enddo
769:                enddo765:                enddo
770:             enddo766:             enddo
771:             enddo767:             enddo
772:          enddo768:          enddo
773:       enddo769:       enddo
811:             !print *, 'coords 708: ', x(1:3*natoms)807:             !print *, 'coords 708: ', x(1:3*natoms)
812:             xdum(:) = x(:)808:             xdum(:) = x(:)
813:             call cart2frac_rb_tri(nrigidbody, xdum, x, H_inverse)809:             call cart2frac_rb_tri(nrigidbody, xdum, x, H_inverse)
814:             !if (ortho) call cart2frac_rb_ortho(xdum, x)810:             !if (ortho) call cart2frac_rb_ortho(xdum, x)
815:             !print *, 'coords 711: ', x(1:3*natoms)811:             !print *, 'coords 711: ', x(1:3*natoms)
816:          endif812:          endif
817:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XR, X)813:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XR, X)
818:          X(:) = XR(:)814:          X(:) = XR(:)
819:       ENDIF815:       ENDIF
820: 816: 
821:       if (boxderivt) call constrain_volume(v_fact, dv_fact, energy1, box_paramsgrad(4:6), gtest)817:       call constrain_volume(v_fact, dv_fact, energy1, box_paramsgrad(4:6), gtest)
822: 818: 
823:       ENERGY = (ENERGY1 + ENERGY2 + ENERGY3)*2625.499D0819:       ENERGY = (ENERGY1 + ENERGY2 + ENERGY3)*2625.499D0
824:       IF (GTEST) G(:) = (G(:) + G3(:))*2625.499D0820:       IF (GTEST) G(:) = (G(:) + G3(:))*2625.499D0
825:       if (gtest) box_paramsgrad(1:6) = box_paramsgrad(1:6)*2625.499D0821:       if (gtest) box_paramsgrad(1:6) = box_paramsgrad(1:6)*2625.499D0
826:       !print *, 'coords benzrigid    : ', x(1:3*natoms)822:       !print *, 'coords benzrigid    : ', x(1:3*natoms)
827:       !print *, 'box params benzrigid: ', box_params(1:6)823:       !print *, 'box params benzrigid: ', box_params(1:6)
828: 824: 
829:       !print *, 'energy: ', energy!1, energy2, energy3825:       !print *, 'energy: ', energy1, energy2, energy3
830:       !print *, 'grad: ', g(1)826:       !print *, 'grad  : '
831:       !!print *, 'grad  : ' 
832:       !do j1 = 1, 2*nrigidbody827:       !do j1 = 1, 2*nrigidbody
833:       !   print *, x(3*j1-2:3*j1)828:       !   print *, g(3*j1-2:3*j1)
834:       !enddo829:       !enddo
835: 830: 
836:       !print *, 'contri: ', energy1, energy2, energy3831:       !print *, 'contri: ', energy1, energy2, energy3
837:       !print *, 'energy: ', energy !, box_params(1:6)832:       !print *, 'energy: ', energy !, box_params(1:6)
838:       !print *, 'g1    : ', g(1:3)833:       !print *, 'g1    : ', g(1:3)
839:       !stop834:       !stop
840: 835: 
841:       END SUBROUTINE BENZGENRIGIDEWALD836:       END SUBROUTINE BENZGENRIGIDEWALD
842: 837: 
843: !     ----------------------------------------------------------------------------------------------838: !     ----------------------------------------------------------------------------------------------


r33131/benzgenrigid_ortho.f90 2017-08-07 15:30:27.133211904 +0100 r33130/benzgenrigid_ortho.f90 2017-08-07 15:30:28.917235637 +0100
  1: ! dj337: Anisotropic potential for polycyclic aromatic hydrocarbons.  1: ! dj337: Anisotropic potential for polycyclic aromatic hydrocarbons.
  2: ! Long-range electrostatic interactions are computed using Ewald summation.  2: ! Long-range electrostatic interactions are computed using Ewald summation.
  3: ! Implemented within the GENRIGID framework.  3: ! Implemented within the GENRIGID framework.
  4:   4: 
  5:       SUBROUTINE BENZGENRIGIDEWALD_ORTHO(X, G, ENERGY, GTEST)  5:       SUBROUTINE BENZGENRIGIDEWALD(X, G, ENERGY, GTEST)
  6:   6: 
  7:       USE COMMONS, ONLY: NATOMS, NCARBON, RBSTLA, RHOCC0, RHOCC10, RHOCC20, &  7:       USE COMMONS, ONLY: NATOMS, NCARBON, RBSTLA, RHOCC0, RHOCC10, RHOCC20, &
  8:      &                   RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, RHOCH10, RHOC20H, &  8:      &                   RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, RHOCH10, RHOC20H, &
  9:      &                   RHOCH20, ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, &  9:      &                   RHOCH20, ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, &
 10:      &                   EWALDREALC, BOX_PARAMS, BOX_PARAMSGRAD 10:      &                   EWALDREALC, BOX_PARAMS, BOX_PARAMSGRAD
 11:  11: 
 12:       ! dj337: PAHA adapted to the genrigid framework 12:       ! dj337: PAHA adapted to the genrigid framework
 13:       USE GENRIGID, ONLY: NRIGIDBODY, ATOMRIGIDCOORDT, TRANSFORMCTORIGID, NSITEPERBODY, & 13:       USE GENRIGID, ONLY: NRIGIDBODY, ATOMRIGIDCOORDT, TRANSFORMCTORIGID, NSITEPERBODY, &
 14:      &                    MAXSITE, SITESRIGIDBODY, TRANSFORMRIGIDTOC, TRANSFORMGRAD 14:      &                    MAXSITE, SITESRIGIDBODY, TRANSFORMRIGIDTOC, TRANSFORMGRAD, INVERSEMATRIX
 15:  15: 
 16:       ! dj337: use Ewald summation to compute electrostatics 16:       ! dj337: use Ewald summation to compute electrostatics
 17:       USE EWALD 17:       USE EWALD
 18:       USE CARTDIST 18:       USE CARTDIST
 19:       USE BOX_DERIVATIVES 19:       USE BOX_DERIVATIVES
 20:  20: 
 21:       IMPLICIT NONE 21:       IMPLICIT NONE
 22:  22: 
 23:       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, IDX
 24:       INTEGER          :: NEWALDREAL(3) 24:       INTEGER          :: NEWALDREAL(3)
 25:       DOUBLE PRECISION :: X(3*NATOMS) 25:       DOUBLE PRECISION :: X(3*NATOMS)
 26:       DOUBLE PRECISION, INTENT(OUT) :: G(3*NATOMS) 26:       DOUBLE PRECISION, INTENT(OUT) :: G(3*NATOMS)
 27:       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), graddum(3*natoms)
 28:       DOUBLE PRECISION, INTENT(OUT) :: ENERGY 28:       DOUBLE PRECISION, INTENT(OUT) :: ENERGY
 29:       DOUBLE PRECISION :: R2, R6, ABSRIJ, DVDR, ENERGY1, ENERGY2, ENERGY3 29:       DOUBLE PRECISION :: R2, R6, ABSRIJ, DVDR, ENERGY1, ENERGY2, ENERGY3, diff, eplus, eminus
 30:       DOUBLE PRECISION :: DMPFCT_SHIFT, EXPFCT_SHIFT, VSHIFT1, VSHIFT2, EWALDREALC2 30:       DOUBLE PRECISION :: DMPFCT_SHIFT, EXPFCT_SHIFT, VSHIFT1, VSHIFT2, EWALDREALC2, RCOMMIN(3), RCOM(3)
 31:       DOUBLE PRECISION :: RI(3), RR(3), RSS(3), NR(3), P(3), EI(3), EJ(3), FRIJ(3), TIJ(3), TJI(3)  31:       DOUBLE PRECISION :: RI(3), RR(3), RSS(3), RSSMIN(3), NR(3), P(3), EI(3), EJ(3), FRIJ(3), TIJ(3), TJI(3) 
 32:       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)
 33:       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)
 34:       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)
 35:       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)
 36:       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 
 37:       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)
  38:       DOUBLE PRECISION :: H(3,3), H_grad(3,3,6), H_inverse(3,3), rrfrac(3), rssfracmin(3), rrcom(3), rcomfrac(3)
  39:       double precision :: rrcomfrac(3), rcomfracmin(3), rssfrac(3), vol, v_fact, dv_fact(3), c(3), s(3)
  40:       double precision :: H_mat(3,3), rmatrix(3), reciplatvec(3,3), reciplatvec_grad(3,3,6)
 38:       DOUBLE PRECISION, PARAMETER :: B = 1.6485D0 41:       DOUBLE PRECISION, PARAMETER :: B = 1.6485D0
  42:       double precision, parameter :: pi = 3.141592654d0
 39:       integer, parameter          :: image_cutoff = 5 43:       integer, parameter          :: image_cutoff = 5
 40:       LOGICAL          :: GTEST 44:       LOGICAL          :: GTEST, keep_angles
 41:  45: 
 42:       !print *, 'hiya' 46:       !print *, 'box_params: ', box_params(1:6)
 43:  47: 
 44:       !if (boxderivt) then 48:       keep_angles = check_angles(box_params(4:6))
 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:  49: 
 53:       !call build_H(H, H_grad, gtest) 50:       if (.not.keep_angles) then
 54:       !call inversematrix(H, H_inverse) 51:          !print *, 'rejecting'
 55:       !call get_volume(vol) 52:          call reject(energy, g)
 56:       !call get_reciplatvec(reciplatvec,reciplatvec_grad, .false.) 53:          return
  54:       endif
 57:  55: 
 58:       newaldreal(:) = floor(ewaldrealc/box_params(1:3) + 0.5d0) 56:       call build_H(H, H_grad, gtest)
 59:       !newaldreal(1) = floor(ewaldrealc*dsqrt(sum(reciplatvec(1,:)**2))/(2.0d0*pi) + 0.5d0) 57:       call inversematrix(H, H_inverse)
 60:       !newaldreal(2) = floor(ewaldrealc*dsqrt(sum(reciplatvec(2,:)**2))/(2.0d0*pi) + 0.5d0) 58:       call get_volume(vol)
 61:       !newaldreal(3) = floor(ewaldrealc*dsqrt(sum(reciplatvec(3,:)**2))/(2.0d0*pi) + 0.5d0) 59:       call get_reciplatvec(reciplatvec,reciplatvec_grad, .false.)
  60: 
  61:       newaldreal(1) = floor(ewaldrealc*dsqrt(sum(reciplatvec(1,:)**2))/(2.0d0*pi) + 0.5d0)
  62:       newaldreal(2) = floor(ewaldrealc*dsqrt(sum(reciplatvec(2,:)**2))/(2.0d0*pi) + 0.5d0)
  63:       newaldreal(3) = floor(ewaldrealc*dsqrt(sum(reciplatvec(3,:)**2))/(2.0d0*pi) + 0.5d0)
 62:       !print *, 'newaldreal: ', newaldreal(:3) 64:       !print *, 'newaldreal: ', newaldreal(:3)
 63:  65: 
 64:       if (boxderivt) then 66:       if (.not. all(newaldreal.le.image_cutoff)) then
 65:          if (.not. all(newaldreal.le.image_cutoff)) then 67:          !print *, 'rejecting'
 66:             print *, 'rejecting2' 68:          call reject(energy, g)
 67:             call reject(energy, g) 69:          return
 68:             return 
 69:          endif 
 70:       endif 70:       endif
 71:  71: 
 72:       ! factorials 72:       ! factorials
 73:       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
 74:       ! initialize energy values 74:       ! initialize energy values
 75:       ! energy1 is due to short-range anisotropic interactions 75:       ! energy1 is due to short-range anisotropic interactions
 76:       ! energy2 is due to damped dispersion 76:       ! energy2 is due to damped dispersion
 77:       ! energy3 is due to long-range electrostatics (computed using Ewald) 77:       ! energy3 is due to long-range electrostatics (computed using Ewald)
 78:       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
 79:  79: 
 85:       !print *, 'coords benzrigid   : ', x(1:25) 85:       !print *, 'coords benzrigid   : ', x(1:25)
 86:  86: 
 87:       ! dj337: check if input coordinates are cartesian 87:       ! dj337: check if input coordinates are cartesian
 88:       ! assumes ATOMRIGIDCOORDT is correct 88:       ! assumes ATOMRIGIDCOORDT is correct
 89:       IF (ATOMRIGIDCOORDT) THEN ! if input is cartesian 89:       IF (ATOMRIGIDCOORDT) THEN ! if input is cartesian
 90:          !print *, 'converting...' 90:          !print *, 'converting...'
 91:          ! convert to rigidbody coordinates 91:          ! convert to rigidbody coordinates
 92:          XR(:) = 0.D0 92:          XR(:) = 0.D0
 93:          CALL TRANSFORMCTORIGID(X, XR) 93:          CALL TRANSFORMCTORIGID(X, XR)
 94:          if (boxderivt) then 94:          if (boxderivt) then
 95:             call frac2cart_rb_ortho(nrigidbody, xdum, xr) 95:             call frac2cart_rb_tri(nrigidbody, xdum, xr, H)
 96:             !if (ortho) call frac2cart_rb_ortho(xdum, xr) 96:             !if (ortho) call frac2cart_rb_ortho(xdum, xr)
 97:             x(:) = xdum(:) 97:             x(:) = xdum(:)
 98:          else 98:          else
 99:             x(:) = xr(:) 99:             x(:) = xr(:)
100:          endif100:          endif
101:       ENDIF101:       ENDIF
102: 102: 
103:       !print *, 'coords benzrigid   : '103:       !print *, 'coords benzrigid   : '
104:       !do j1 = 1, 2*nrigidbody104:       !do j1 = 1, 2*nrigidbody
105:       !   print *, x(3*j1-2:3*j1)105:       !   print *, x(3*j1-2:3*j1)
106:       !enddo106:       !enddo
107: 107: 
108:       !call build_H(H, H_grad, gtest)108:       !call build_H(H, H_grad, gtest)
109:       !call inversematrix(H, H_inverse)109:       !call inversematrix(H, H_inverse)
110:       !call get_volume(vol)110:       !call get_volume(vol)
111:       !call get_reciplatvec(reciplatvec,reciplatvec_grad, .false.)111:       !call get_reciplatvec(reciplatvec,reciplatvec_grad, .false.)
112:       !if (boxderivt) then112:       ! compute v factor
113:       !   ! compute v factor113:       c(:) = dcos(box_params(4:6))
114:       !   c(:) = dcos(box_params(4:6))114:       s(:) = dsin(box_params(4:6))
115:       !   s(:) = dsin(box_params(4:6))115:       v_fact = dsqrt(1.0d0 - c(1)**2-c(2)**2-c(3)**2 + 2.0d0*c(1)*c(2)*c(3))
116:       !   v_fact = dsqrt(1.0d0 - c(1)**2-c(2)**2-c(3)**2 + 2.0d0*c(1)*c(2)*c(3))116:       dv_fact(1) = s(1)*(c(1) - c(2)*c(3))/v_fact
117:       !   dv_fact(1) = s(1)*(c(1) - c(2)*c(3))/v_fact117:       dv_fact(2) = s(2)*(c(2) - c(1)*c(3))/v_fact
118:       !   dv_fact(2) = s(2)*(c(2) - c(1)*c(3))/v_fact118:       dv_fact(3) = s(3)*(c(3) - c(1)*c(2))/v_fact
119:       !   dv_fact(3) = s(3)*(c(3) - c(1)*c(2))/v_fact 
120:       !endif 
121: 119: 
122:       EWALDREALC2 = EWALDREALC**2120:       EWALDREALC2 = EWALDREALC**2
123: 121: 
124:       !print *, 'reciplatvec: ', reciplatvec(:3,:3)122:       !print *, 'reciplatvec: ', reciplatvec(:3,:3)
125: 123: 
126:       ! OFFSET is number of CoM coords (3*NRIGIDBODY)124:       ! OFFSET is number of CoM coords (3*NRIGIDBODY)
127:       OFFSET     = 3*NRIGIDBODY125:       OFFSET     = 3*NRIGIDBODY
128: 126: 
129:       ! Computing Cartesian coordinates for the system.  127:       ! Computing Cartesian coordinates for the system.  
130:       DO J1 = 1, NRIGIDBODY128:       DO J1 = 1, NRIGIDBODY
198:                J6 = OFFSET + J4196:                J6 = OFFSET + J4
199: 197: 
200:                ! loop over sites in the rigid body J2198:                ! loop over sites in the rigid body J2
201:                DO J = 1, NSITEPERBODY(J2)199:                DO J = 1, NSITEPERBODY(J2)
202: 200: 
203:                   ! J8 is index for site J201:                   ! J8 is index for site J
204:                   J8     = MAXSITE*(J2-1) + J202:                   J8     = MAXSITE*(J2-1) + J
205:                   ! EJ is Z-axis for site J203:                   ! EJ is Z-axis for site J
206:                   EJ(:)  = E(J8,:)204:                   EJ(:)  = E(J8,:)
207:                   rr(:) = r(j7,:) - r(j8,:)205:                   rr(:) = r(j7,:) - r(j8,:)
208:                   ! minimum image convention 
209:                   rssmin(1) = rr(1) - box_params(1)*anint(rr(1)/box_params(1)) 
210:                   rssmin(2) = rr(2) - box_params(2)*anint(rr(2)/box_params(2)) 
211:                   rssmin(3) = rr(3) - box_params(3)*anint(rr(3)/box_params(3)) 
212:                   ! convert to fractional coordinates206:                   ! convert to fractional coordinates
213:                   !rrfrac(:) = matmul(H_inverse, rr(:))207:                   rrfrac(:) = matmul(H_inverse, rr(:))
214:                   !! minimum image convention208:                   ! minimum image convention
215:                   !rssfracmin(1) = rrfrac(1) - anint(rrfrac(1))209:                   rssfracmin(1) = rrfrac(1) - anint(rrfrac(1))
216:                   !rssfracmin(2) = rrfrac(2) - anint(rrfrac(2))210:                   rssfracmin(2) = rrfrac(2) - anint(rrfrac(2))
217:                   !rssfracmin(3) = rrfrac(3) - anint(rrfrac(3))211:                   rssfracmin(3) = rrfrac(3) - anint(rrfrac(3))
218: 212: 
219:                   if (gtest.and.boxderivt) then213:                   if (gtest.and.boxderivt) then
220:                      ! get center of mass separation vector214:                      ! get center of mass separation vector
221:                      rrcom(:) = x(j3-2:j3) - x(j4-2:j4)215:                      rrcom(:) = x(j3-2:j3) - x(j4-2:j4)
222:                      ! minimum image convention 
223:                      rrcommin(1) = rrcom(1) - box_params(1)*anint(rr(1)/box_params(1)) 
224:                      rrcommin(2) = rrcom(2) - box_params(2)*anint(rr(2)/box_params(2)) 
225:                      rrcommin(3) = rrcom(3) - box_params(3)*anint(rr(3)/box_params(3)) 
226:                      ! convert to fractional coordinates216:                      ! convert to fractional coordinates
227:                      !rrcomfrac(:) = matmul(H_inverse, rrcom(:))217:                      rrcomfrac(:) = matmul(H_inverse, rrcom(:))
228:                      !! minimum image convention218:                      ! minimum image convention
229:                      !rcomfracmin(1) = rrcomfrac(1) - anint(rrfrac(1))219:                      rcomfracmin(1) = rrcomfrac(1) - anint(rrfrac(1))
230:                      !rcomfracmin(2) = rrcomfrac(2) - anint(rrfrac(2))220:                      rcomfracmin(2) = rrcomfrac(2) - anint(rrfrac(2))
231:                      !rcomfracmin(3) = rrcomfrac(3) - anint(rrfrac(3))221:                      rcomfracmin(3) = rrcomfrac(3) - anint(rrfrac(3))
232:                   endif222:                   endif
233: 223: 
234:                   do l = -newaldreal(1), newaldreal(1)224:                   do l = -newaldreal(1), newaldreal(1)
235:                   rss(1) = rssmin(1) + box_params(1)*l225:                   rssfrac(1) = rssfracmin(1) + l
236: 226: 
237:                      do m = -newaldreal(2), newaldreal(2)227:                      do m = -newaldreal(2), newaldreal(2)
238:                      rss(2) = rssmin(2) + box_params(2)*m228:                      rssfrac(2) = rssfracmin(2) + m
239: 229: 
240:                         do n = -newaldreal(3), newaldreal(3)230:                         do n = -newaldreal(3), newaldreal(3)
241:                         rss(3) = rssmin(3) + box_params(3)*n231:                         rssfrac(3) = rssfracmin(3) + n
242: 232: 
243:                         ! convert to absolute coordinates233:                         ! convert to absolute coordinates
244:                         !rss(:) = matmul(H, rssfrac(:))234:                         rss(:) = matmul(H, rssfrac(:))
245: 235: 
246:                         ! get COM vector236:                         ! get COM vector
247:                         if (gtest.and.boxderivt) then237:                         if (gtest.and.boxderivt) then
248:                            rcom(1) = rrcommin(1) + box_params(1)*l238:                            rcomfrac(1) = rcomfracmin(1) + l
249:                            rcom(2) = rrcommin(2) + box_params(2)*m239:                            rcomfrac(2) = rcomfracmin(2) + m
250:                            rcom(3) = rrcommin(3) + box_params(3)*n240:                            rcomfrac(3) = rcomfracmin(3) + n
 241:                            ! convert to absolute coordinates
 242:                            rcom(:) = matmul(H, rcomfrac(:))
251:                         endif243:                         endif
252:                      244:                      
253:                         R2     = DOT_PRODUCT(RSS(:),RSS(:))245:                         R2     = DOT_PRODUCT(RSS(:),RSS(:))
254:                         !print *, 'r2: ', r2246:                         !print *, 'r2: ', r2
255:                         ! check if distance within cutoff247:                         ! check if distance within cutoff
256:                         IF (R2 < EWALDREALC2) THEN248:                         IF (R2 < EWALDREALC2) THEN
257:                            !print *, j7, j8249:                            !print *, j7, j8
258:                            !print *, 'r   : ', rss(:3)250:                            !print *, 'r   : ', rss(:3)
259:                            !print *, 'rmin: ', rssmin(1:3)251:                            !print *, 'rmin: ', rssmin(1:3)
260:                            !print *, 'rcom: ', rcom(1:3)252:                            !print *, 'rcom: ', rcom(1:3)
466:                               !g(j5-2:j5) = g(j5-2:j5) + tij(:)458:                               !g(j5-2:j5) = g(j5-2:j5) + tij(:)
467:                               ! total gradient wrt AA coords for rigid body J2459:                               ! total gradient wrt AA coords for rigid body J2
468:                               G(J6-2:J6) = G(J6-2:J6) + DVDR*DRIJDPJ(:) + TJI(:)460:                               G(J6-2:J6) = G(J6-2:J6) + DVDR*DRIJDPJ(:) + TJI(:)
469:                               !g(j6-2:j6) = g(j6-2:j6) + tji(:)461:                               !g(j6-2:j6) = g(j6-2:j6) + tji(:)
470: 462: 
471:                               ! TODO: this is if orientation of rigid body depends on cell parameters463:                               ! 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)464:                               !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(:)))465:                               !rmatrix(:) = matmul(H_inverse, vol*(rss(:)-rcom(:)))
474:                               !box_paramsgrad(4) = box_paramsgrad(4) + dvdr*dot_product(rss(1:3), matmul(H_mat, rmatrix(:)))466:                               !box_paramsgrad(4) = box_paramsgrad(4) + dvdr*dot_product(rss(1:3), matmul(H_mat, rmatrix(:)))
475: 467: 
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(:)))468:                               !box_paramsgrad(1) = box_paramsgrad(1) + dvdr*dot_product(rss(1:3), matmul(H_mat/vol, rmatrix(:)))
478:                               !if (boxderivt) then469: 
479:                               !   do idx = 1, 6470:                               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(:)))471:                                  box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product((dvdr*rss(1:3) + frij(1:3)), matmul(H_grad(:,:,idx), rcomfrac(:)))
481:                               !   enddo472:                               enddo
482:                               !endif 
483: 473: 
484:                            ENDIF474:                            ENDIF
485:       475:       
486:                         ENDIF476:                         ENDIF
487: 477: 
488:                      enddo478:                      enddo
489:                   enddo479:                   enddo
490:                enddo480:                enddo
491: 481: 
492:                ENDDO482:                ENDDO
514:             do j = 1, nsiteperbody(j1)504:             do j = 1, nsiteperbody(j1)
515:                j8 = maxsite*(j1-1) + j505:                j8 = maxsite*(j1-1) + j
516:                ej(:) = e(j8,:)506:                ej(:) = e(j8,:)
517: 507: 
518:                !print *, 'new!!', j1, i, j508:                !print *, 'new!!', j1, i, j
519:                rr(:) = r(j7,:) - r(j8,:)509:                rr(:) = r(j7,:) - r(j8,:)
520:                !print *, j7, r(j7,:3)510:                !print *, j7, r(j7,:3)
521:                !print *, j8, r(j8,:3)511:                !print *, j8, r(j8,:3)
522:                !print *, 'rr: ', rr(:3)512:                !print *, 'rr: ', rr(:3)
523:                ! convert to fractional513:                ! convert to fractional
524:                !rrfrac(:) = matmul(H_inverse, rr(:))514:                rrfrac(:) = matmul(H_inverse, rr(:))
525:                !print *, 'rrfrac: ', rrfrac(:3)515:                !print *, 'rrfrac: ', rrfrac(:3)
526: 516: 
527:                ! loop over boxes517:                ! loop over boxes
528:                do l = -newaldreal(1), newaldreal(1)518:                do l = -newaldreal(1), newaldreal(1)
529:                   do m = -newaldreal(2), newaldreal(2)519:                   do m = -newaldreal(2), newaldreal(2)
530:                      do n = -newaldreal(3), newaldreal(3)520:                      do n = -newaldreal(3), newaldreal(3)
531: 521: 
532:                      if (.not.(l.eq.0.and.m.eq.0.and.n.eq.0)) then522:                      if (.not.(l.eq.0.and.m.eq.0.and.n.eq.0)) then
533: 523: 
534:                         rss(1) = rr(1) + box_params(1)*l524:                         rssfrac(1) = rrfrac(1) + l
535:                         rss(2) = rr(2) + box_params(2)*m525:                         rssfrac(2) = rrfrac(2) + m
536:                         rss(3) = rr(3) + box_params(3)*n526:                         rssfrac(3) = rrfrac(3) + n
537:                         !print *, l, m, n527:                         !print *, l, m, n
538:                         !print *, 'rssfrac: ', rssfrac(:3)528:                         !print *, 'rssfrac: ', rssfrac(:3)
539:                         !rss(:) = matmul(H, rssfrac(:))529:                         rss(:) = matmul(H, rssfrac(:))
540: 530: 
541:                         if (gtest.and.boxderivt) then531:                         if (gtest.and.boxderivt) then
542:                            rcom(1) = box_params(1)*l532:                            rcomfrac(1) = l
543:                            rcom(2) = box_params(2)*m533:                            rcomfrac(2) = m
544:                            rcom(3) = box_params(3)*n534:                            rcomfrac(3) = n
 535:                            rcom(:) = matmul(H, rcomfrac(:))
545:                         endif536:                         endif
546: 537: 
547:                         r2 = dot_product(rss(:), rss(:))538:                         r2 = dot_product(rss(:), rss(:))
548:                         !print *, 'rssmin2: ', rss(1:3)539:                         !print *, 'rssmin2: ', rss(1:3)
549:                         if (r2 < ewaldrealc2) then540:                         if (r2 < ewaldrealc2) then
550: 541: 
551:                         !print *, 'rssmin3: ', rss(1:3)542:                         !print *, 'rssmin3: ', rss(1:3)
552:                         !print *, 'r2    : ', r2543:                         !print *, 'r2    : ', r2
553:                         !print *, 'rcom  : ', rcom(1:3)544:                         !print *, 'rcom  : ', rcom(1:3)
554:                         absrij = dsqrt(r2)545:                         absrij = dsqrt(r2)
564:                         DMPFCT = 1.D0555:                         DMPFCT = 1.D0
565:                         DMPFCT_SHIFT = 1.D0556:                         DMPFCT_SHIFT = 1.D0
566:                         ! initialize sum for the derivative of damping function557:                         ! initialize sum for the derivative of damping function
567:                         DDMPDR = B558:                         DDMPDR = B
568: 559: 
569:                         ! calculate sums560:                         ! calculate sums
570:                         DO K = 1, 6561:                         DO K = 1, 6
571: 562: 
572:                            DMPFCT = DMPFCT + (B*ABSRIJ)**K/FLOAT(FCT(K))563:                            DMPFCT = DMPFCT + (B*ABSRIJ)**K/FLOAT(FCT(K))
573:                            DMPFCT_SHIFT = DMPFCT_SHIFT + (B*EWALDREALC)**K/FLOAT(FCT(K))564:                            DMPFCT_SHIFT = DMPFCT_SHIFT + (B*EWALDREALC)**K/FLOAT(FCT(K))
574:                            IF (K > 1) DDMPDR = DDMPDR + (B**K)*(ABSRIJ)**(K-1)/FLOAT(FCT(K-1))565:                            !IF (K > 1) DDMPDR = DDMPDR + (B**K)*(ABSRIJ)**(K-1)/FLOAT(FCT(K-1))
575: 566: 
576:                         END DO567:                         END DO
577: 568: 
578:                         EXPFCT = DEXP(-B*ABSRIJ)569:                         EXPFCT = DEXP(-B*ABSRIJ)
579:                         EXPFCT_SHIFT = DEXP(-B*EWALDREALC)570:                         EXPFCT_SHIFT = DEXP(-B*EWALDREALC)
580:                         ! DDMPDR is derivative of damping function with factor 1/Rab571:                         ! DDMPDR is derivative of damping function with factor 1/Rab
581:                         DDMPDR = (B*EXPFCT*DMPFCT - EXPFCT*DDMPDR)/ABSRIJ572:                         DDMPDR = (B*EXPFCT*DMPFCT - EXPFCT*DDMPDR)/ABSRIJ
582:                         ! DMPFCT is damping function573:                         ! DMPFCT is damping function
583:                         DMPFCT = 1.D0 - EXPFCT*DMPFCT574:                         DMPFCT = 1.D0 - EXPFCT*DMPFCT
584:                         ! DMPFCT_SHIFT is vertical shift for damping function575:                         ! DMPFCT_SHIFT is vertical shift for damping function
755:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPJ(:) + TJI(:)746:                            G(J5-2:J5) = G(J5-2:J5) + DVDR*DRIJDPJ(:) + TJI(:)
756:                            !g(j5-2:j5) = g(j5-2:j5) + tji(:)747:                            !g(j5-2:j5) = g(j5-2:j5) + tji(:)
757:                            !print *, 'dispersion:'748:                            !print *, 'dispersion:'
758:                            !print *, dvdr*drijdpi(:3)749:                            !print *, dvdr*drijdpi(:3)
759:                            !print *, dvdr*drijdpj(:3)750:                            !print *, dvdr*drijdpj(:3)
760:                            !print *, 'energy2: ', energy2751:                            !print *, 'energy2: ', energy2
761:                            !print *, 'anisotropic:'752:                            !print *, 'anisotropic:'
762:                            !print *, tij(:3)753:                            !print *, tij(:3)
763:                            !print *, tji(:3)754:                            !print *, tji(:3)
764: 755: 
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)756:                            do idx = 1, 6
766:                            !if (boxderivt) then757:                                  box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product((dvdr*rss(1:3) + frij(1:3)), matmul(H_grad(:,:,idx), rcomfrac(:)))
767:                            !   do idx = 1, 6758:                                  !box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product(dvdr*rss(1:3), matmul(H_grad(:,:,idx), rcomfrac(:)))
768:                            !         box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product((dvdr*rss(1:3) + frij(1:3)), matmul(H_grad(:,:,idx), rcomfrac(:)))759:                            enddo
769:                            !         !box_paramsgrad(idx) = box_paramsgrad(idx) + dot_product(dvdr*rss(1:3), matmul(H_grad(:,:,idx), rcomfrac(:))) 
770:                            !   enddo 
771:                            !endif 
772: 760: 
773:                         ENDIF761:                         ENDIF
774:                         endif762:                         endif
775:                     endif763:                     endif
776:                   enddo764:                   enddo
777:                enddo765:                enddo
778:             enddo766:             enddo
779:             enddo767:             enddo
780:          enddo768:          enddo
781:       enddo769:       enddo
782: 770: 
783:       ! convert to cartesian coordinates771:       ! convert to cartesian coordinates
784:       XC(:) = 0.D0772:       XC(:) = 0.D0
785:       if (boxderivt) then773:       if (boxderivt) then
786:          xdum(:) = x(:)774:          xdum(:) = x(:)
787:          call cart2frac_rb_ortho(nrigidbody, xdum, x)775:          call cart2frac_rb_tri(nrigidbody, xdum, x, H_inverse)
788:       endif776:       endif
789:       CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XC, X)777:       CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XC, X)
790:       ! restore cartesian rigid body coordinates778:       ! restore cartesian rigid body coordinates
791:       if (boxderivt) x(:) = xdum(:)779:       if (boxderivt) x(:) = xdum(:)
792: 780: 
793: !      !!! ENERGY3 and G3 are energy and gradient due to electrostatics781: !      !!! ENERGY3 and G3 are energy and gradient due to electrostatics
794: !      !!! computed using Ewald summation782: !      !!! computed using Ewald summation
795:       CALL EWALDSUM(1, XC, G3C, ENERGY3, GTEST)783:       CALL EWALDSUM(1, XC, G3C, ENERGY3, GTEST)
796: !784: !
797: !      !!! convert Ewald contribution of gradient to rigidbody coordinates785: !      !!! convert Ewald contribution of gradient to rigidbody coordinates
811:       !print *, 'energy: ', energy2799:       !print *, 'energy: ', energy2
812: 800: 
813:       ! dj337: if input was cartesian, convert back to cartesian801:       ! dj337: if input was cartesian, convert back to cartesian
814:       ! assumes ATOMRIGIDCOORDT is correct802:       ! assumes ATOMRIGIDCOORDT is correct
815:       IF (ATOMRIGIDCOORDT) THEN803:       IF (ATOMRIGIDCOORDT) THEN
816: 804: 
817:          ! convert to cartesian coordinates805:          ! convert to cartesian coordinates
818:          if (boxderivt) then806:          if (boxderivt) then
819:             !print *, 'coords 708: ', x(1:3*natoms)807:             !print *, 'coords 708: ', x(1:3*natoms)
820:             xdum(:) = x(:)808:             xdum(:) = x(:)
821:             call cart2frac_rb_ortho(nrigidbody, xdum, x)809:             call cart2frac_rb_tri(nrigidbody, xdum, x, H_inverse)
822:             !if (ortho) call cart2frac_rb_ortho(xdum, x)810:             !if (ortho) call cart2frac_rb_ortho(xdum, x)
823:             !print *, 'coords 711: ', x(1:3*natoms)811:             !print *, 'coords 711: ', x(1:3*natoms)
824:          endif812:          endif
825:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XR, X)813:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XR, X)
826:          X(:) = XR(:)814:          X(:) = XR(:)
827:       ENDIF815:       ENDIF
828: 816: 
829:       !if (boxderivt) call constrain_volume(v_fact, dv_fact, energy1, box_paramsgrad(4:6), gtest)817:       call constrain_volume(v_fact, dv_fact, energy1, box_paramsgrad(4:6), gtest)
830: 818: 
831:       ENERGY = (ENERGY1 + ENERGY2 + ENERGY3)*2625.499D0819:       ENERGY = (ENERGY1 + ENERGY2 + ENERGY3)*2625.499D0
832:       IF (GTEST) G(:) = (G(:) + G3(:))*2625.499D0820:       IF (GTEST) G(:) = (G(:) + G3(:))*2625.499D0
833:       if (gtest) box_paramsgrad(1:3) = box_paramsgrad(1:3)*2625.499D0821:       if (gtest) box_paramsgrad(1:6) = box_paramsgrad(1:6)*2625.499D0
834:       !print *, 'coords benzrigid    : ', x(1:3*natoms)822:       !print *, 'coords benzrigid    : ', x(1:3*natoms)
835:       !print *, 'box params benzrigid: ', box_params(1:6)823:       !print *, 'box params benzrigid: ', box_params(1:6)
836: 824: 
837:       !print *, 'energy: ', energy!1, energy2, energy3825:       !print *, 'energy: ', energy1, energy2, energy3
838:       !print *, 'grad: ', g(1)826:       !print *, 'grad  : '
839:       !!print *, 'grad  : ' 
840:       !do j1 = 1, 2*nrigidbody827:       !do j1 = 1, 2*nrigidbody
841:       !   print *, x(3*j1-2:3*j1)828:       !   print *, g(3*j1-2:3*j1)
842:       !enddo829:       !enddo
843: 830: 
844:       !print *, 'contri: ', energy1, energy2, energy3831:       !print *, 'contri: ', energy1, energy2, energy3
845:       !print *, 'energy: ', energy !, box_params(1:6)832:       !print *, 'energy: ', energy !, box_params(1:6)
846:       !print *, 'g1    : ', g(1:3)833:       !print *, 'g1    : ', g(1:3)
847:       !stop834:       !stop
848: 835: 
849:       END SUBROUTINE BENZGENRIGIDEWALD_ORTHO836:       END SUBROUTINE BENZGENRIGIDEWALD
850: 837: 
851: !     ----------------------------------------------------------------------------------------------838: !     ----------------------------------------------------------------------------------------------
 839: !
 840: !      SUBROUTINE DEFPAHARIGID()
 841: !
 842: !      USE COMMONS, ONLY: RHOCC0, RHOCC10, RHOCC20,  RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, RHOCH10, RHOC20H, RHOCH20, &
 843: !                         ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, CCKJ
 844: !
 845: !      IMPLICIT NONE
 846: ! 
 847: !      ALPHACC = 1.861500D0
 848: !      ALPHAHH = 1.431200D0
 849: !      ALPHACH = 1.775600D0
 850: !
 851: !      DC6CC    = 30.469D0
 852: !      DC6HH    = 5.359D0
 853: !      DC6CH    = 12.840D0
 854: !
 855: !      RHOCC0  = 5.814700D0
 856: !      RHOCC10 = 0.021700D0
 857: !      RHOCC20 =-0.220800D0
 858: !
 859: !      RHOHH0  = 4.486200D0
 860: !      RHOHH10 =-0.271800D0
 861: !      RHOHH20 = 0.0D0
 862: !
 863: !      RHOCH0  = 5.150500D0
 864: !      RHOC10H = 0.021700D0
 865: !      RHOCH10 =-0.271800D0
 866: !      RHOC20H =-0.220800D0
 867: !      RHOCH20 = 0.0D0
 868: !
 869: !      KKJ     = 1.D-03
 870: !      CCKJ    = 1.D0   !1389.354848D0
 871: !
 872: !      END SUBROUTINE DEFPAHARIGID
 873: !
 874: !!     ----------------------------------------------------------------------------------------------
 875: 
 876:       SUBROUTINE DEFBENZENERIGIDEWALD()
 877: 
 878:       USE COMMONS, ONLY: NRBSITES, SITE, RBSTLA, STCHRG 
 879:       USE GENRIGID, ONLY: NRIGIDBODY
 880: 
 881:       IMPLICIT NONE
 882:  
 883:       INTEGER :: J1
 884: 
 885: !!     C6H6
 886: 
 887: !!     D6h reference geometry: C-C: 1.397 angstrom; C-H: 1.087 angstrom
 888: !      ! dj337: note that all done in atomic units
 889: 
 890:       SITE(1,:)  = (/ 2.63923430843701,   0.00000000000000,   0.00000000000000/)
 891:       SITE(2,:)  = (/ 1.31961715421850,  -2.28564395764590,   0.00000000000000/)
 892:       SITE(3,:)  = (/-1.31961715421850,  -2.28564395764590,   0.00000000000000/)
 893:       SITE(4,:)  = (/-2.63923430843701,   0.00000000000000,   0.00000000000000/)
 894:       SITE(5,:)  = (/-1.31961715421850,   2.28564395764590,   0.00000000000000/)
 895:       SITE(6,:)  = (/ 1.31961715421850,   2.28564395764590,   0.00000000000000/)
 896:       SITE(7,:)  = (/ 4.69338981379532,   0.00000000000000,   0.00000000000000/)
 897:       SITE(8,:)  = (/ 2.34669490689766,  -4.06459480860986,   0.00000000000000/)
 898:       SITE(9,:)  = (/-2.34669490689766,  -4.06459480860986,   0.00000000000000/)
 899:       SITE(10,:) = (/-4.69338981379532,   0.00000000000000,   0.00000000000000/)
 900:       SITE(11,:) = (/-2.34669490689766,   4.06459480860986,   0.00000000000000/)
 901:       SITE(12,:) = (/ 2.34669490689766,   4.06459480860986,   0.00000000000000/)
 902: 
 903:       RBSTLA(1,:)  = SITE(7,:)  - SITE(1,:)                 ! Z FROM C1 TO H1
 904:       RBSTLA(2,:)  = SITE(8,:)  - SITE(2,:)                 ! Z FROM C2 TO H2
 905:       RBSTLA(3,:)  = SITE(9,:)  - SITE(3,:)                 ! Z FROM C3 TO H3
 906:       RBSTLA(4,:)  = SITE(10,:) - SITE(4,:)                 ! Z FROM C4 TO H4
 907:       RBSTLA(5,:)  = SITE(11,:) - SITE(5,:)                 ! Z FROM C5 TO H5
 908:       RBSTLA(6,:)  = SITE(12,:) - SITE(6,:)                 ! Z FROM C6 TO H6
 909:       RBSTLA(7,:)  = SITE(7,:)  - SITE(1,:)                 ! Z FROM C1 TO H1!
 910:       RBSTLA(8,:)  = SITE(8,:)  - SITE(2,:)                 ! Z FROM C2 TO H2!
 911:       RBSTLA(9,:)  = SITE(9,:) -  SITE(3,:)                 ! Z FROM C3 TO H3!
 912:       RBSTLA(10,:) = SITE(10,:) - SITE(4,:)                 ! Z FROM C4 TO H4!
 913:       RBSTLA(11,:) = SITE(11,:) - SITE(5,:)                 ! Z FROM C5 TO H5!
 914:       RBSTLA(12,:) = SITE(12,:) - SITE(6,:)                 ! Z FROM C6 TO H6!
 915: 
 916:       DO J1 = 1, NRBSITES
 917:  
 918:          RBSTLA(J1,:)   = RBSTLA(J1,:)/DSQRT(DOT_PRODUCT(RBSTLA(J1,:),RBSTLA(J1,:)))
 919: 
 920:       ENDDO
 921: 
 922:       DO J1 = 1, NRIGIDBODY
 923:          STCHRG(12*(J1-1)+1:12*(J1-1)+6) = -0.11114D0
 924:          STCHRG(12*(J1-1)+7:12*(J1-1)+12) = 0.11114D0
 925:       ENDDO
 926: 
 927:       END SUBROUTINE DEFBENZENERIGIDEWALD


r33131/box_derivatives.f90 2017-08-07 15:30:27.353214831 +0100 r33130/box_derivatives.f90 2017-08-07 15:30:29.133238511 +0100
125:  125:  
126: ! -----------------------------------------------------------------------------------126: ! -----------------------------------------------------------------------------------
127: 127: 
128: ! CONVERTS GRADIENT WRT ATOMIC POSITIONS FROM ABSOLUTE TO FRACTIONAL COORDINATES.128: ! CONVERTS GRADIENT WRT ATOMIC POSITIONS FROM ABSOLUTE TO FRACTIONAL COORDINATES.
129: ! Wrapper subroutine that calls appropriate subroutine depending on whether atomistic129: ! Wrapper subroutine that calls appropriate subroutine depending on whether atomistic
130: ! or rigidbody coordinates, orthorhombic or triclinic unit cell.130: ! or rigidbody coordinates, orthorhombic or triclinic unit cell.
131: 131: 
132:       subroutine boxderiv(grad, gradfrac, H)132:       subroutine boxderiv(grad, gradfrac, H)
133: 133: 
134:       use commons, only: ortho134:       use commons, only: ortho
135:       use genrigid, only: rigidinit, atomrigidcoordt135:       use genrigid, only: rigidinit, atomrigidcoordt, inversematrix
136: 136: 
137:       implicit none137:       implicit none
138: 138: 
139:       double precision, intent(in)           :: grad(3*natoms)139:       double precision, intent(in)           :: grad(3*natoms)
140:       double precision, intent(out)          :: gradfrac(3*natoms)140:       double precision, intent(out)          :: gradfrac(3*natoms)
141:       double precision, intent(in), optional :: H(3,3)141:       double precision, intent(in), optional :: H(3,3)
142:       double precision                       :: H_mat(3,3), H_grad(3,3,6)142:       double precision                       :: H_mat(3,3), H_grad(3,3,6)
143: 143: 
144:       ! orthorhombic cell144:       ! orthorhombic cell
145:       if (ortho) then145:       if (ortho) then


r33131/checkd.f90 2017-08-07 15:30:27.577217811 +0100 r33130/checkd.f90 2017-08-07 15:30:29.353241438 +0100
  1:       SUBROUTINE CHECKD(X)  1:       SUBROUTINE CHECKD(X)
  2:   2: 
  3:       USE COMMONS, ONLY: NATOMS, COMPRESST, PERCOLATET, CHECKDID, GTHOMSONT, &  3:       USE COMMONS, ONLY: NATOMS, COMPRESST, PERCOLATET, CHECKDID, GTHOMSONT, &
  4:                          BOXDERIVT, ORTHO, BOX_PARAMS, BOX_PARAMSGRAD  4:                          BOXDERIVT, ORTHO, BOX_PARAMS, BOX_PARAMSGRAD
  5:       USE GENRIGID, ONLY: RIGIDINIT, ATOMRIGIDCOORDT, DEGFREEDOMS, TRANSFORMCTORIGID  5:       USE GENRIGID, ONLY: RIGIDINIT, ATOMRIGIDCOORDT, DEGFREEDOMS, TRANSFORMCTORIGID
  6:   6: 
  7:       USE MODHESS  7:       USE MODHESS
  8:       IMPLICIT NONE  8:       IMPLICIT NONE
  9:   9: 
 10:       INTEGER          :: IVRNO, IVRNO1, IVRNO2, J1, J3, dof, doff 10:       INTEGER          :: IVRNO, IVRNO1, IVRNO2, J1, J3
 11:       DOUBLE PRECISION :: X(3*NATOMS), G(3*NATOMS), ENERGY, FM, FP, DFA, DFN, TMPCOORDS(3*NATOMS) 11:       DOUBLE PRECISION :: X(3*NATOMS), G(3*NATOMS), ENERGY, FM, FP, DFA, DFN, TMPCOORDS(3*NATOMS)
 12:       double precision :: box_paramsold(6) 12:       double precision :: box_paramsold(6)
 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:       !print *, 'checkd coords: '
 18:       if (rigidinit) then 18:       !do j1 = 1, natoms
 19:          dof = degfreedoms 19:       !   print *, x(3*j1-2:3*j1)
 20:       else 20:       !enddo
 21:          dof = 3*natoms 
 22:       endif 
 23:  21: 
 24:       print *, 'DELX: ', DELX 22:       print *, 'DELX: ', DELX
 25:  23: 
 26: ! jwrm2> Turning compression on, if required 24: ! jwrm2> Turning compression on, if required
 27:       IF (COMPRESST .OR. PERCOLATET) COMPON = .TRUE. 25:       IF (COMPRESST .OR. PERCOLATET) COMPON = .TRUE.
 28:  26: 
 29: ! jwrm2> Converting GTHOMSON coordinates to polars 27: ! jwrm2> Converting GTHOMSON coordinates to polars
 30:       IF (GTHOMSONT) THEN 28:       IF (GTHOMSONT) THEN
 31:         CALL GTHOMSONCTOANG(X(1:3*NATOMS), TMPCOORDS(1:3*NATOMS), NATOMS) 29:         CALL GTHOMSONCTOANG(X(1:3*NATOMS), TMPCOORDS(1:3*NATOMS), NATOMS)
 32:         X(1:3*NATOMS) = TMPCOORDS(1:3*NATOMS) 30:         X(1:3*NATOMS) = TMPCOORDS(1:3*NATOMS)
 36:  34: 
 37:       IF (CHECKDID == 0) THEN 35:       IF (CHECKDID == 0) THEN
 38:          GTEST = .FALSE. 36:          GTEST = .FALSE.
 39:          CALL POTENTIAL (X, G, ENERGY, GTEST, STEST) 37:          CALL POTENTIAL (X, G, ENERGY, GTEST, STEST)
 40:          WRITE(*, *) 'Energy  = ', ENERGY 38:          WRITE(*, *) 'Energy  = ', ENERGY
 41:  39: 
 42:       ELSEIF (CHECKDID == 1) THEN 40:       ELSEIF (CHECKDID == 1) THEN
 43:  41: 
 44: !     Checks gradients 42: !     Checks gradients
 45:  43: 
 46:       ! check derivatives wrt atomic positions 44:       ! dj337: if box derivatives are computed
 47:       DO IVRNO = 1, DOF 45:       IF (BOXDERIVT) THEN
 48:          WRITE(*, *) IVRNO 46:          !IF (ORTHO) THEN
 49:  47:             ! check derivatives wrt atomic positions
 50:          ! dj337: make sure coordinates rigid body 48:             DO IVRNO = 1, DEGFREEDOMS!3*NATOMS
 51:          if (rigidinit.and.atomrigidcoordt) then 49:                WRITE(*, *) IVRNO
 52:             call transformctorigid(x, tmpcoords) 50: 
 53:             x(1:degfreedoms) = tmpcoords(1:degfreedoms) 51:                if (rigidinit.and.atomrigidcoordt) then
 54:             x(degfreedoms+1:3*natoms) = 0.0d0 52:                   !print *, 'oh heyyyy'
 55:             atomrigidcoordt = .false. 53:                   call transformctorigid(x, tmpcoords)
 56:          endif 54:                   x(1:degfreedoms) = tmpcoords(1:degfreedoms)
 57:  55:                   x(degfreedoms+1:3*natoms) = 0.0d0
 58:          !print *, 'x: ', x(ivrno) 56:                   atomrigidcoordt = .false.
 59:  57:                endif
 60:          GTEST = .FALSE. 58: 
 61:          X(IVRNO) = X(IVRNO) - DELX 59:                !print *, x(1:3*natoms)
 62:          !print *, 'xplus1: ', x(ivrno) 60:                GTEST = .FALSE.
 63:          CALL POTENTIAL(X, G, FM, GTEST, STEST) 61:                X(IVRNO) = X(IVRNO) - DELX
 64:          !print *, 'xplus2: ', x(ivrno) 62:                !write(*, *) 'X minus  = ', x(ivrno)
 65:  63:                CALL POTENTIAL(X, G, FM, GTEST, STEST)
 66:          if (rigidinit.and.atomrigidcoordt) then 64:                !write(*, *) 'X minus2 = ', x(ivrno)
 67:             call transformctorigid(x, tmpcoords) 65:                !WRITE(*, *) 'Energy minus = ', FM
 68:             x(1:degfreedoms) = tmpcoords(1:degfreedoms) 66: 
 69:             x(degfreedoms+1:3*natoms) = 0.0d0 67:                if (rigidinit.and.atomrigidcoordt) then
 70:             atomrigidcoordt = .false. 68:                   call transformctorigid(x, tmpcoords)
 71:          endif 69:                   x(1:degfreedoms) = tmpcoords(1:degfreedoms)
 72:  70:                   x(degfreedoms+1:3*natoms) = 0.0d0
 73:          X(IVRNO) = X(IVRNO) + 2.D0*DELX 71:                   atomrigidcoordt = .false.
 74:          !print *, 'xminus1: ', x(ivrno) 72:                endif
 75:          CALL POTENTIAL(X, G, FP, GTEST, STEST) 73: 
 76:          !print *, 'xminus2: ', x(ivrno) 74:                X(IVRNO) = X(IVRNO) + 2.D0*DELX
  75:                !write(*, *) 'X plus   = ', x(ivrno)
  76:                CALL POTENTIAL(X, G, FP, GTEST, STEST)
  77:                !WRITE(*, *) 'Energy plus  = ', FP
  78:                !write(*, *) 'X plus2  = ', x(ivrno)
 77:       79:      
 78:          if (rigidinit.and.atomrigidcoordt) then 80:                if (rigidinit.and.atomrigidcoordt) then
 79:             call transformctorigid(x, tmpcoords) 81:                   call transformctorigid(x, tmpcoords)
 80:             x(1:degfreedoms) = tmpcoords(1:degfreedoms) 82:                   x(1:degfreedoms) = tmpcoords(1:degfreedoms)
 81:             x(degfreedoms+1:3*natoms) = 0.0d0 83:                   x(degfreedoms+1:3*natoms) = 0.0d0
 82:             atomrigidcoordt = .false.  84:                   atomrigidcoordt = .false. 
 83:          endif 85:                endif
 84:   86:  
 85:          GTEST = .TRUE. 87:                GTEST = .TRUE.
 86:          X(IVRNO) = X(IVRNO) - DELX 88:                X(IVRNO) = X(IVRNO) - DELX
 87:          CALL POTENTIAL(X, G, ENERGY, GTEST, STEST) 89:                !write(*, *) 'X       = ', x(ivrno)
 88:          DFN = (FP - FM) / (2.D0*DELX) 90:                CALL POTENTIAL(X, G, ENERGY, GTEST, STEST)
 89:          IF (ABS(DFN).LT.1.0D-10) DFN = 0.D0 91:                !write(*, *) 'X2      = ', x(ivrno)
 90:          DFA = G(IVRNO) 92:                DFN = (FP - FM) / (2.D0*DELX)
 91:  93:                IF (ABS(DFN).LT.1.0D-10) DFN = 0.D0
 92:          WRITE(*, *) 'Gradient numerical  = ', DFN 94:                DFA = G(IVRNO)
 93:          WRITE(*, *) 'Gradient analytical = ', DFA 95: 
 94:  96:                WRITE(*, *) 'Gradient numerical  = ', DFN
 95:          IF (ABS(DFN - DFA) > ERRLIM) WRITE(*, *) IVRNO, DFN, DFA, ABS(DFN-DFA) 97:                WRITE(*, *) 'Gradient analytical = ', DFA
 96:       ENDDO 98: 
 97:  99:                IF (ABS(DFN - DFA) > ERRLIM) WRITE(*, *) IVRNO, DFN, DFA, ABS(DFN-DFA)
 98:       ! dj337: check lattice derivatives100:                !stop
 99:       if (boxderivt) then101:             ENDDO
100:          doff = 3102: 
101:          if (.not.ortho) doff = 6103:             print *, 'New system: '
102:          DO IVRNO = 1, doff104:             ! check lattice derivatives
103: 105:             DO IVRNO = 1,6
104:             if (rigidinit.and.atomrigidcoordt) then106: 
105:                call transformctorigid(x, tmpcoords)107:                if (rigidinit.and.atomrigidcoordt) then
106:                x(1:degfreedoms) = tmpcoords(1:degfreedoms)108:                   call transformctorigid(x, tmpcoords)
107:                atomrigidcoordt = .false.109:                   x(1:degfreedoms) = tmpcoords(1:degfreedoms)
108:             endif110:                   atomrigidcoordt = .false.
109: 111:                endif
110:             WRITE(*, *) 'Box parameter ', IVRNO112: 
111: 113:                WRITE(*, *) 'Box parameter ', IVRNO
112:             GTEST = .FALSE.114: 
113:             BOX_PARAMS(IVRNO) = BOX_PARAMS(IVRNO) - DELX115:                GTEST = .FALSE.
114:             CALL POTENTIAL(X, G, FM, GTEST, STEST)116:                !box_paramsold(:) = box_params(:)
 117:                !print *, 'box_paramsold: ', box_paramsold(1:6)
 118:                BOX_PARAMS(IVRNO) = BOX_PARAMS(IVRNO) - DELX
 119:                !print *, 'atomrigidcoordt: ', atomrigidcoordt
 120:                !call rotate_bodies(box_paramsold, x)
 121:                !write(*, *) 'Box minus  = ', box_params(ivrno)
 122:                CALL POTENTIAL(X, G, FM, GTEST, STEST)
 123:                !WRITE(*, *) 'Box minus2 = ', BOX_PARAMS(IVRNO)
 124:                WRITE(*, *) 'Energy minus = ', FM
 125: 
 126:                !box_paramsold(:) = box_params(:)
 127:                !print *, 'box_paramsold: ', box_paramsold(1:6)
 128:                BOX_PARAMS(IVRNO) = BOX_PARAMS(IVRNO) + 2.D0*DELX
 129:                !print *, 'atomrigidcoordt: ', atomrigidcoordt
 130:                !call rotate_bodies(box_paramsold, x)
 131:                !write(*, *), 'Box plus  = ', box_params(ivrno)
 132:                CALL POTENTIAL(X, G, FP, GTEST, STEST)
 133:                !WRITE(*, *) 'Box plus2  = ', BOX_PARAMS(IVRNO)
 134:                WRITE(*, *) 'Energy plus  = ', FP
 135: 
 136:                GTEST = .TRUE.
 137:                !box_paramsold(:) = box_params(:)
 138:                !print *, 'box_paramsold: ', box_paramsold(1:6)
 139:                BOX_PARAMS(IVRNO) = BOX_PARAMS(IVRNO) - DELX
 140:                !print *, 'atomrigidcoordt: ', atomrigidcoordt
 141:                !call rotate_bodies(box_paramsold, x)
 142:                !write(*, *), 'Box       = ', box_params(ivrno) 
 143:                CALL POTENTIAL(X, G, ENERGY, GTEST, STEST)
 144:                !write(*, *), 'Box2      = ', box_params(ivrno)
 145:                DFN = (FP - FM) / (2.D0*DELX)
 146:                IF (ABS(DFN).LT.1.0D-10) DFN = 0.D0
 147:                DFA = BOX_PARAMSGRAD(IVRNO)
 148: 
 149:                WRITE(*, *) 'Box gradient numerical  = ', DFN
 150:                WRITE(*, *) 'Box gradient analytical = ', DFA
 151: 
 152:                IF (ABS(DFN - DFA) > ERRLIM) WRITE(*, *) IVRNO, DFN, DFA, ABS(DFN-DFA)
 153: 
 154:             ENDDO
 155:          !ENDIF
 156: 
 157:       ELSE ! if no box derivatives
 158: 
 159:          DO IVRNO = 1, 3*NATOMS!DEGFREEDOMS
 160: 
 161:             WRITE(*, *) IVRNO
 162: 
 163:             GTEST    = .FALSE.
 164:             X(IVRNO) = X(IVRNO) - DELX
 165:             CALL POTENTIAL (X, G, FM, GTEST, STEST)
115:             WRITE(*, *) 'Energy minus = ', FM166:             WRITE(*, *) 'Energy minus = ', FM
116: 167: 
117:             BOX_PARAMS(IVRNO) = BOX_PARAMS(IVRNO) + 2.D0*DELX168:             X(IVRNO) = X(IVRNO) + 2.D0*DELX
118:             CALL POTENTIAL(X, G, FP, GTEST, STEST)169:             CALL POTENTIAL (X, G,  FP, GTEST, STEST)
119:             WRITE(*, *) 'Energy plus  = ', FP170:             WRITE(*, *) 'Energy plus  = ', FP
120: 171: 
121:             GTEST = .TRUE.172:             GTEST = .TRUE.
122:             BOX_PARAMS(IVRNO) = BOX_PARAMS(IVRNO) - DELX173:             X(IVRNO) = X(IVRNO) - DELX
123:             CALL POTENTIAL(X, G, ENERGY, GTEST, STEST)174:             CALL POTENTIAL (X, G, ENERGY, GTEST, STEST)
124:             DFN = (FP - FM) / (2.D0*DELX)175:             DFN = (FP - FM) / (2.D0*DELX)
125:             IF (ABS(DFN).LT.1.0D-10) DFN = 0.D0176:             IF (ABS(DFN) .LT. 1.0D-10) DFN = 0.D0
126:             DFA = BOX_PARAMSGRAD(IVRNO)177:             DFA = G(IVRNO)
127: 178: 
128:             WRITE(*, *) 'Box gradient numerical  = ', DFN179:             WRITE(*, *) 'Gradient numerical  = ', DFN
129:             WRITE(*, *) 'Box gradient analytical = ', DFA180:             WRITE(*, *) 'Gradient analytical = ', DFA
130: 181: 
131:             IF (ABS(DFN - DFA) > ERRLIM) WRITE(*, *) IVRNO, DFN, DFA, ABS(DFN-DFA)182:             IF (ABS(DFN - DFA) > ERRLIM) WRITE(*, *) IVRNO, DFN, DFA, ABS(DFN-DFA)
132: 183: 
133:          ENDDO184:          ENDDO
134:       endif185: 
 186:       ENDIF
135: 187: 
136:       ELSE IF (CHECKDID == 2) THEN188:       ELSE IF (CHECKDID == 2) THEN
137: 189: 
138:          IF (.NOT. ALLOCATED(HESS)) ALLOCATE(HESS(3*NATOMS,3*NATOMS))190:          IF (.NOT. ALLOCATED(HESS)) ALLOCATE(HESS(3*NATOMS,3*NATOMS))
139: 191: 
140:          DO IVRNO1 = 1, 3*NATOMS192:          DO IVRNO1 = 1, 3*NATOMS
141:             DO IVRNO2 = 1, 3*NATOMS193:             DO IVRNO2 = 1, 3*NATOMS
142:                WRITE(*,*) IVRNO1, IVRNO2194:                WRITE(*,*) IVRNO1, IVRNO2
143:                X(IVRNO1) = X(IVRNO1) - DELX195:                X(IVRNO1) = X(IVRNO1) - DELX
144:                CALL POTENTIAL (X,G,ENERGY,.TRUE.,.FALSE.)196:                CALL POTENTIAL (X,G,ENERGY,.TRUE.,.FALSE.)


r33131/ewald.f90 2017-08-07 15:30:27.801220791 +0100 r33130/ewald.f90 2017-08-07 15:30:29.581244472 +0100
511:       use commons, only: rerhoarray, imrhoarray 511:       use commons, only: rerhoarray, imrhoarray 
512:       use cartdist, only: get_reciplatvec512:       use cartdist, only: get_reciplatvec
513: 513: 
514:       implicit none514:       implicit none
515: 515: 
516:       integer                      :: j1, j3, l, m, n, dims(3)516:       integer                      :: j1, j3, l, m, n, dims(3)
517:       integer, intent(in)          :: newaldrecip(3)517:       integer, intent(in)          :: newaldrecip(3)
518:       double precision, intent(in) :: x(3*natoms)518:       double precision, intent(in) :: x(3*natoms)
519:       double precision             :: k(3), r(3), reciplatvec(3,3), reciplatvec_grad(3,3,6)519:       double precision             :: k(3), r(3), reciplatvec(3,3), reciplatvec_grad(3,3,6)
520:       double precision             :: q1, k2, kdotr, rerho, imrho, ewaldrecipc2520:       double precision             :: q1, k2, kdotr, rerho, imrho, ewaldrecipc2
 521:       double precision, parameter  :: pi = 3.141592654D0
521: 522: 
522:       ! reciprocal-space cutoff523:       ! reciprocal-space cutoff
523:       ewaldrecipc2 = ewaldrecipc**2524:       ewaldrecipc2 = ewaldrecipc**2
524: 525: 
525:       ! make sure allocated arrays for structure factors are the correct size526:       ! make sure allocated arrays for structure factors are the correct size
526:       dims(:) = 2*newaldrecip(1:3)+1 527:       dims(:) = 2*newaldrecip(1:3)+1 
527:       if (.not.allocated(rerhoarray)) allocate(rerhoarray(dims(1), dims(2), dims(3)))528:       if (.not.allocated(rerhoarray)) allocate(rerhoarray(dims(1), dims(2), dims(3)))
528:       if (.not.allocated(imrhoarray)) allocate(imrhoarray(dims(1), dims(2), dims(3)))529:       if (.not.allocated(imrhoarray)) allocate(imrhoarray(dims(1), dims(2), dims(3)))
529: 530: 
530:       if (.not.(size(rerhoarray,1).eq.dims(1).and.size(rerhoarray,2).eq.dims(2).and.size(rerhoarray,3).eq.dims(3))) then531:       if (.not.(size(rerhoarray,1).eq.dims(1).and.size(rerhoarray,2).eq.dims(2).and.size(rerhoarray,3).eq.dims(3))) then


r33131/keywords.f 2017-08-07 15:30:28.025223770 +0100 r33130/keywords.f 2017-08-07 15:30:29.813247555 +0100
5421:             BOXLZ=BOXLX5421:             BOXLZ=BOXLX
5422:             BOX3D(3) = BOX3D(1)5422:             BOX3D(3) = BOX3D(1)
5423:          ENDIF5423:          ENDIF
5424:          PI = 4.D0*DATAN(1.D0)5424:          PI = 4.D0*DATAN(1.D0)
5425:          PALPHA=PI/2.D05425:          PALPHA=PI/2.D0
5426:          PBETA=PALPHA5426:          PBETA=PALPHA
5427:          PGAMMA=PALPHA5427:          PGAMMA=PALPHA
5428:          IF (NITEMS.GT.4) THEN5428:          IF (NITEMS.GT.4) THEN
5429:             ORTHO = .FALSE.5429:             ORTHO = .FALSE.
5430:             CALL READF(XX)5430:             CALL READF(XX)
5431:             PGAMMA=XX5431:             PGAMMA=PI*XX/180.D0
5432:          ENDIF5432:          ENDIF
5433:          IF (NITEMS.GT.6) THEN5433:          IF (NITEMS.GT.6) THEN
5434:             PALPHA=PGAMMA5434:             PALPHA=PGAMMA
5435:             CALL READF(XX)5435:             CALL READF(XX)
5436:             PBETA=XX5436:             PBETA=PI*XX/180.D0
5437:             CALL READF(XX)5437:             CALL READF(XX)
5438:             PGAMMA=XX5438:             PGAMMA=PI*XX/180.D0
5439:          ENDIF5439:          ENDIF
5440: 5440: 
5441:          ! dj337: also put into box_params vector5441:          ! dj337: also put into box_params vector
5442:          box_params(1) = boxlx5442:          box_params(1) = boxlx
5443:          box_params(2) = boxly5443:          box_params(2) = boxly
5444:          box_params(3) = boxlz5444:          box_params(3) = boxlz
5445:          box_params(4) = palpha5445:          box_params(4) = palpha
5446:          box_params(5) = pbeta5446:          box_params(5) = pbeta
5447:          box_params(6) = pgamma5447:          box_params(6) = pgamma
5448: ! lower triangular lattice matrix5448: ! lower triangular lattice matrix


r33131/main.F 2017-08-07 15:30:28.249226750 +0100 r33130/main.F 2017-08-07 15:30:30.033250483 +0100
374: 374: 
375:       IF (GAUSST) THEN375:       IF (GAUSST) THEN
376:          ALLOCATE(GAUSSKK(3*NATOMSALLOC,GMODES),GAUSSEE(GMODES),GKSMALL(3*NATOMSALLOC))376:          ALLOCATE(GAUSSKK(3*NATOMSALLOC,GMODES),GAUSSEE(GMODES),GKSMALL(3*NATOMSALLOC))
377:          CALL KEGEN ! INITIAL SETUP 377:          CALL KEGEN ! INITIAL SETUP 
378:          DO J1=1,GMODES378:          DO J1=1,GMODES
379:              PRINT *,J1,GAUSSEE(J1)379:              PRINT *,J1,GAUSSEE(J1)
380:          ENDDO380:          ENDDO
381:       ENDIF381:       ENDIF
382: 382: 
383:       QMINP(1:NSAVE,1:3*NATOMSALLOC)=0.0D0 ! to prevent reading from uninitialised memory383:       QMINP(1:NSAVE,1:3*NATOMSALLOC)=0.0D0 ! to prevent reading from uninitialised memory
384:       if (boxderivt) boxq(1:nsave,1:6) = 0.0d0384:       boxq(1:nsave,1:6) = 0.0d0
385:       QMINT(1:NSAVE,1:NATOMSALLOC)=1 ! to prevent reading from uninitialised memory385:       QMINT(1:NSAVE,1:NATOMSALLOC)=1 ! to prevent reading from uninitialised memory
386:       QMINNATOMS(1:NSAVE)=NATOMS ! to prevent reading from uninitialised memory386:       QMINNATOMS(1:NSAVE)=NATOMS ! to prevent reading from uninitialised memory
387:       AVOIDNATOMS(1:MAXSAVE)=NATOMS ! as above387:       AVOIDNATOMS(1:MAXSAVE)=NATOMS ! as above
388:       COORDSO(1:3*NATOMSALLOC,1:NPAR)=0.0D0 ! to prevent reading from uninitialised memory388:       COORDSO(1:3*NATOMSALLOC,1:NPAR)=0.0D0 ! to prevent reading from uninitialised memory
389:       FF(1:NSAVE)=0 ! to prevent reading from uninitialised memorY389:       FF(1:NSAVE)=0 ! to prevent reading from uninitialised memorY
390:       VATO(1:NATOMSALLOC,1:NPAR)=0.0D0 ! to prevent reading from uninitialised memory390:       VATO(1:NATOMSALLOC,1:NPAR)=0.0D0 ! to prevent reading from uninitialised memory
391:       ALLOCATE(ESAVE(NTAB,NPAR),XINSAVE(NTAB,NPAR))391:       ALLOCATE(ESAVE(NTAB,NPAR),XINSAVE(NTAB,NPAR))
392:       ALLOCATE(VEC(NVEC))392:       ALLOCATE(VEC(NVEC))
393: 393: 
394: !     IF (SYMMETRIZE.AND.(.NOT.CENT)) THEN394: !     IF (SYMMETRIZE.AND.(.NOT.CENT)) THEN


r33131/potential.f90 2017-08-07 15:30:28.473229730 +0100 r33130/potential.f90 2017-08-07 15:30:30.253253412 +0100
353: !         EMINUS=EDUM1+EDUM2353: !         EMINUS=EDUM1+EDUM2
354: !354: !
355: !         X(J1)=X(J1)+DIFF355: !         X(J1)=X(J1)+DIFF
356: !         GRAD(J1)=(EPLUS-EMINUS)/(2.0D0*DIFF)356: !         GRAD(J1)=(EPLUS-EMINUS)/(2.0D0*DIFF)
357: !         WRITE(*,'(I5, 4F20.10)') J1, GRAD(J1), EPLUS, EMINUS, ABS(EPLUS-EMINUS)357: !         WRITE(*,'(I5, 4F20.10)') J1, GRAD(J1), EPLUS, EMINUS, ABS(EPLUS-EMINUS)
358: !     END DO358: !     END DO
359: !359: !
360: 360: 
361: ! dj337: pahagenrigid for benzene with Ewald361: ! dj337: pahagenrigid for benzene with Ewald
362:    ELSE IF (BENZRIGIDEWALDT) THEN362:    ELSE IF (BENZRIGIDEWALDT) THEN
363:       if (ortho) then363:       CALL BENZGENRIGIDEWALD(X, GRAD, EREAL, GRADT)
364:          call benzgenrigidewald_ortho(x, grad, ereal, gradt) 
365:       else 
366:          CALL BENZGENRIGIDEWALD(X, GRAD, EREAL, GRADT) 
367:       endif 
368:       if (rigidinit .and. (.not.atomrigidcoordt)) then364:       if (rigidinit .and. (.not.atomrigidcoordt)) then
369:          xrigidcoords(1:degfreedoms) = x(1:degfreedoms)365:          xrigidcoords(1:degfreedoms) = x(1:degfreedoms)
370:          xrigidgrad(1:degfreedoms) = grad(1:degfreedoms)366:          xrigidgrad(1:degfreedoms) = grad(1:degfreedoms)
371:          call transformrigidtoc(1, nrigidbody, xcoords, xrigidcoords)367:          call transformrigidtoc(1, nrigidbody, xcoords, xrigidcoords)
372:       endif368:       endif
373: ! check analytical and numerical gradients369: ! check analytical and numerical gradients
374: !      DIFF=1.0D-6370: !      DIFF=1.0D-6
375: !      WRITE(*, *) 'analytic and numerical gradients:'371: !      WRITE(*, *) 'analytic and numerical gradients:'
376: !      print *, 'diff: ', diff372: !      print *, 'diff: ', diff
377: !      IF (ATOMRIGIDCOORDT) THEN ! if input is cartesian373: !      IF (ATOMRIGIDCOORDT) THEN ! if input is cartesian


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0