From b08f823f78b13ba05149e7f3275306de581de1dd Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Mon, 25 Nov 2019 12:00:10 +0100 Subject: [PATCH 1/1] working gradient for cations --- source/unres/data/energy_data.F90 | 3 +- source/unres/energy.F90 | 69 ++++++++++++++++++++----------------- source/unres/io_config.F90 | 11 ++++++ 3 files changed, 51 insertions(+), 32 deletions(-) diff --git a/source/unres/data/energy_data.F90 b/source/unres/data/energy_data.F90 index 8f88483..55353dc 100644 --- a/source/unres/data/energy_data.F90 +++ b/source/unres/data/energy_data.F90 @@ -409,7 +409,8 @@ 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 diff --git a/source/unres/energy.F90 b/source/unres/energy.F90 index 91ca96b..44e994d 100644 --- a/source/unres/energy.F90 +++ b/source/unres/energy.F90 @@ -11935,13 +11935,13 @@ ! 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 @@ -11957,8 +11957,8 @@ !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 @@ -22644,7 +22644,7 @@ !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, & @@ -22658,7 +22658,7 @@ 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 @@ -22670,7 +22670,7 @@ 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) @@ -22679,6 +22679,8 @@ 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) @@ -22818,11 +22820,11 @@ 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) @@ -22880,8 +22882,14 @@ 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 @@ -22975,11 +22983,11 @@ 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 @@ -23045,17 +23053,17 @@ 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 @@ -23138,7 +23146,6 @@ !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. @@ -26350,7 +26357,7 @@ 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& @@ -26360,13 +26367,13 @@ - 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)& @@ -26374,7 +26381,7 @@ - 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) & @@ -26906,15 +26913,15 @@ 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 @@ -27241,25 +27248,25 @@ + 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) diff --git a/source/unres/io_config.F90 b/source/unres/io_config.F90 index 5f8c061..2906fb5 100644 --- a/source/unres/io_config.F90 +++ b/source/unres/io_config.F90 @@ -3287,6 +3287,17 @@ ! 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 -- 1.7.9.5