working gradient for cations
authorAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Mon, 25 Nov 2019 11:00:10 +0000 (12:00 +0100)
committerAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Mon, 25 Nov 2019 11:00:10 +0000 (12:00 +0100)
source/unres/data/energy_data.F90
source/unres/energy.F90
source/unres/io_config.F90

index 8f88483..55353dc 100644 (file)
         real(kind=8),dimension(:,:),allocatable :: alphapolcat,&
            epsheadcat,sig0headcat,sigiso1cat,sigiso2cat,rborncat,&
            sigmap1cat,sigmap2cat,chiscat,wquadcat,chippcat,&
-           epsintabcat,debaykapcat,chicat,sigmacat, nstatecat, epscat
+           epsintabcat,debaykapcat,chicat,sigmacat, nstatecat, epscat,&
+           aa_aq_cat,bb_aq_cat
 
         real(kind=8),dimension(:,:,:),allocatable :: alphasurcat,&
            alphisocat,wqdipcat,dtailcat,wstatecat
index 91ca96b..44e994d 100644 (file)
 !       print *,i,j,gg_lipi(3),gg_lipj(3),sss_ele_cut
 !      write (iout,*) "gg",(gg(k),k=1,3)
       do k=1,3
-        gvdwx(k,i)=gvdwx(k,i)-gg(k) +gg_lipi(k)&
+        gradpepcatx(k,i)=gradpepcatx(k,i)-gg(k) &
                   +(eom12*(dc_norm(k,j)-om12*dc_norm(k,nres+i)) &
                   +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
 
-        gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)&
-                  +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,j)) &
-                  +eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv   
+!        gradpepcatx(k,j)=gradpepcatx(k,j)+gg(k) &
+!                  +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,j)) &
+!                  +eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv   
 
 !        write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
 !                 +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
 !grad        enddo
 !grad      enddo
       do l=1,3
-        gvdwc(l,i)=gvdwc(l,i)-gg(l)
-        gvdwc(l,j)=gvdwc(l,j)+gg(l)
+        gradpepcat(l,i)=gradpepcat(l,i)-gg(l)
+        gradpepcat(l,j)=gradpepcat(l,j)+gg(l)
       enddo
       end subroutine sc_grad_cat
 
 !c------------------------------------------------------------------------------
 #endif
       subroutine ecatcat(ecationcation)
-        integer :: i,j,itmp,xshift,yshift,zshift,subchap,k
+        integer :: i,j,itmp,xshift,yshift,zshift,subchap,k,itypi,itypj
         real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,&
         r7,r4,ecationcation,k0,rcal
         real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
         epscalc=0.05
         r06 = rcat0**6
         r012 = r06**2
-        k0 = 332.0*(2.0*2.0)/80.0
+!        k0 = 332.0*(2.0*2.0)/80.0
         itmp=0
         
         do i=1,4
         xi=c(1,i)
         yi=c(2,i)
         zi=c(3,i)
-         
+          itypi=itype(i,5)
           xi=mod(xi,boxxsize)
           if (xi.lt.0) xi=xi+boxxsize
           yi=mod(yi,boxysize)
           if (zi.lt.0) zi=zi+boxzsize
 
           do j=i+1,itmp+nres_molec(5)
+          itypj=itype(j,5)
+          k0 = 332.0*(ichargecat(itypi)*ichargecat(itypj))/80.0
 !           print *,i,j,'catcat'
            xj=c(1,j)
            yj=c(2,j)
          do j=itmp+1,itmp+nres_molec(5)
 
 ! Calculate SC interaction energy.
-            itypj=iabs(itype(j,1))
+            itypj=iabs(itype(j,5))
             if ((itypj.eq.ntyp1)) cycle
              CALL elgrad_init_cat(eheadtail,Egb,Ecl,Elj,Equad,Epol)
 
-            dscj_inv=vbld_inv(j+nres)
+            dscj_inv=0.0
            xj=c(1,j)
            yj=c(2,j)
            zj=c(3,j)
           chi1 = chicat(itypi,itypj)
           chis1 = chiscat(itypi,itypj)
           chip1 = chippcat(itypi,itypj)
+!          chi1=0.0d0
+!          chis1=0.0d0
+!          chip1=0.0d0
+          chi2=0.0
+          chip2=0.0
+          chis2=0.0
 !          chis2 = chis(itypj,itypi)
-!          chis12 = chis1 * chis2
+          chis12 = chis1 * chis2
           sig1 = sigmap1cat(itypi,itypj)
 !          sig2 = sigmap2(itypi,itypj)
 ! alpha factors from Fcav/Gcav
           sigder = -sig * sigsq
           rij_shift = 1.0D0 / rij_shift
           fac       = rij_shift**expon
-          c1        = fac  * fac * aa_aq(itypi,itypj)
+          c1        = fac  * fac * aa_aq_cat(itypi,itypj)
 !          print *,"ADAM",aa_aq(itypi,itypj)
 
 !          c1        = 0.0d0
-          c2        = fac  * bb_aq(itypi,itypj)
+          c2        = fac  * bb_aq_cat(itypi,itypj)
 !          c2        = 0.0d0
           evdwij    = eps1 * eps2rt * eps3rt * ( c1 + c2 )
           eps2der   = eps3rt * evdwij
        facd2 = dtailcat(2,itypi,itypj) * vbld_inv(j+nres)
        DO k = 1, 3
         pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres))
-        gvdwx(k,i) = gvdwx(k,i) &
+        gradpepcatx(k,i) = gradpepcatx(k,i) &
                   - (( dFdR + gg(k) ) * pom)
         pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres))
 !        gvdwx(k,j) = gvdwx(k,j)   &
 !                  + (( dFdR + gg(k) ) * pom)
-        gvdwc(k,i) = gvdwc(k,i)  &
+        gradpepcat(k,i) = gradpepcat(k,i)  &
                   - (( dFdR + gg(k) ) * ertail(k))
-        gvdwc(k,j) = gvdwc(k,j) &
+        gradpepcat(k,j) = gradpepcat(k,j) &
                   + (( dFdR + gg(k) ) * ertail(k))
         gg(k) = 0.0d0
-
+       ENDDO
 !c! Compute head-head and head-tail energies for each state
           isel = iabs(Qi) + iabs(Qj)
           IF (isel.eq.0) THEN
 !c!-------------------------------------------------------------------
 !c! NAPISY KONCOWE
          END DO   ! j
-        END DO    ! iint
        END DO     ! i
 !c      write (iout,*) "Number of loop steps in EGB:",ind
 !c      energy_dec=.false.
         facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j)))
 
         pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
-        gvdwx(k,i) = gvdwx(k,i) &
+        gradpepcatx(k,i) = gradpepcatx(k,i) &
                   - dGCLdR * pom&
                   - dGGBdR * pom&
                   - dGCVdR * pom&
                   - dGLJdR * pom
 
         pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
-        gvdwx(k,j) = gvdwx(k,j)+ dGCLdR * pom&
+        gradpepcatx(k,j) = gradpepcatx(k,j)+ dGCLdR * pom&
                    + dGGBdR * pom+ dGCVdR * pom&
                   + dPOLdR1 * (erhead_tail(k,1)&
       -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j)))&
                   + dPOLdR2 * condor + dGLJdR * pom
 
-        gvdwc(k,i) = gvdwc(k,i)  &
+        gradpepcat(k,i) = gradpepcat(k,i)  &
                   - dGCLdR * erhead(k)&
                   - dGGBdR * erhead(k)&
                   - dGCVdR * erhead(k)&
                   - dPOLdR2 * erhead_tail(k,2)&
                   - dGLJdR * erhead(k)
 
-        gvdwc(k,j) = gvdwc(k,j)         &
+        gradpepcat(k,j) = gradpepcat(k,j)         &
                   + dGCLdR * erhead(k) &
                   + dGGBdR * erhead(k) &
                   + dGCVdR * erhead(k) &
         condor = (erhead_tail(k,2) &
        + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j)))
 
-        gvdwx(k,i) = gvdwx(k,i) &
+        gradpepcatx(k,i) = gradpepcatx(k,i) &
                    - dPOLdR2 * (erhead_tail(k,2) &
        -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))
-        gvdwx(k,j) = gvdwx(k,j)   &
+        gradpepcatx(k,j) = gradpepcatx(k,j)   &
                    + dPOLdR2 * condor
 
-        gvdwc(k,i) = gvdwc(k,i) &
+        gradpepcat(k,i) = gradpepcat(k,i) &
                    - dPOLdR2 * erhead_tail(k,2)
-        gvdwc(k,j) = gvdwc(k,j) &
+        gradpepcat(k,j) = gradpepcat(k,j) &
                    + dPOLdR2 * erhead_tail(k,2)
 
        END DO
        + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j)))
 
         pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
-        gvdwx(k,i) = gvdwx(k,i) &
+        gradpepcatx(k,i) = gradpepcatx(k,i) &
                   - dGCLdR * pom &
                   - dPOLdR2 * (erhead_tail(k,2) &
        -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) &
                   - dGLJdR * pom
 
         pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
-        gvdwx(k,j) = gvdwx(k,j) &
+        gradpepcatx(k,j) = gradpepcatx(k,j) &
                   + dGCLdR * pom &
                   + dPOLdR2 * condor &
                   + dGLJdR * pom
 
 
-        gvdwc(k,i) = gvdwc(k,i) &
+        gradpepcat(k,i) = gradpepcat(k,i) &
                   - dGCLdR * erhead(k) &
                   - dPOLdR2 * erhead_tail(k,2) &
                   - dGLJdR * erhead(k)
 
-        gvdwc(k,j) = gvdwc(k,j) &
+        gradpepcat(k,j) = gradpepcat(k,j) &
                   + dGCLdR * erhead(k) &
                   + dPOLdR2 * erhead_tail(k,2) &
                   + dGLJdR * erhead(k)
index 5f8c061..2906fb5 100644 (file)
 !       print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i) 
        END DO
       END DO
+      allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
+      do i=1,ntyp
+        do j=1,ntyp_molec(5)
+          epsij=epscat(i,j)
+          rrij=sigmacat(i,j)
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_aq_cat(i,j)=epsij*rrij*rrij
+          bb_aq_cat(i,j)=-sigeps*epsij*rrij
+         enddo
+       enddo 
       endif