X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.F90;h=1e0ed70e1b040726426ca38f6023c0e7415f250a;hb=c6f305ff4a0b1fdc111572dc9692d866b9ffbb15;hp=5c9aacf99464db97144f92f75ab04ce394a5046f;hpb=f8ee39d07f69eabbbf4945b63c93308c85fe5a96;p=unres4.git diff --git a/source/unres/energy.F90 b/source/unres/energy.F90 index 5c9aacf..1e0ed70 100644 --- a/source/unres/energy.F90 +++ b/source/unres/energy.F90 @@ -850,8 +850,11 @@ else call ecatcat(ecationcation) endif + if (oldion.gt.0) then call ecat_prot(ecation_prot) - call ecats_prot_amber(ecations_prot_amber) + else + call ecats_prot_amber(ecation_prot) + endif if (nres_molec(2).gt.0) then call eprot_sc_base(escbase) call epep_sc_base(epepbase) @@ -934,7 +937,7 @@ energia(47)=epepbase energia(48)=escpho energia(49)=epeppho - energia(50)=ecations_prot_amber +! energia(50)=ecations_prot_amber call sum_energy(energia,.true.) if (dyn_ss) call dyn_set_nss ! print *," Processor",myrank," left SUM_ENERGY" @@ -1061,7 +1064,7 @@ epepbase=energia(47) escpho=energia(48) epeppho=energia(49) - ecations_prot_amber=energia(50) +! ecations_prot_amber=energia(50) ! energia(41)=ecation_prot ! energia(42)=ecationcation @@ -1080,7 +1083,7 @@ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl& +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl& +wcatprot*ecation_prot+wcatcat*ecationcation+wscbase*escbase& - +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ecations_prot_amber + +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc & @@ -1094,7 +1097,7 @@ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl& +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl& +wcatprot*ecation_prot+wcatcat*ecationcation+wscbase*escbase& - +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ecations_prot_amber + +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho #endif energia(0)=etot ! detecting NaNQ @@ -1281,7 +1284,7 @@ epepbase=energia(47) escpho=energia(48) epeppho=energia(49) - ecations_prot_amber=energia(50) +! ecations_prot_amber=energia(50) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,& estr,wbond,ebe,wang,& @@ -1297,7 +1300,7 @@ etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,& ecorr3_nucl,wcorr3_nucl,ecation_prot,wcatprot,ecationcation,wcatcat, & escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,& - ecations_prot_amber,etot + etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ & @@ -1360,7 +1363,7 @@ etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,& ecorr3_nucl,wcorr3_nucl,ecation_prot,wcatprot,ecationcation,wcatcat, & escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,& - ecations_prot_amber,etot + etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ & @@ -11893,15 +11896,8 @@ end subroutine sc_grad subroutine sc_grad_cat -! implicit real*8 (a-h,o-z) use calc_data -! include 'DIMENSIONS' -! include 'COMMON.CHAIN' -! include 'COMMON.DERIV' -! include 'COMMON.CALC' -! include 'COMMON.IOUNITS' real(kind=8), dimension(3) :: dcosom1,dcosom2 -! print *,"wchodze" eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 & +dCAVdOM1+ dGCLdOM1+ dPOLdOM1 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 & @@ -11915,11 +11911,6 @@ ! eom2=0.0d0 ! eom12=evdwij*eps1_om12 ! end diagnostics -! write (iout,*) "eps2der",eps2der," eps3der",eps3der,& -! " sigder",sigder -! write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12 -! write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12 -!C print *,sss_ele_cut,'in sc_grad' do k=1,3 dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k)) @@ -11932,13 +11923,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 @@ -11948,17 +11939,44 @@ ! ! Calculate the components of the gradient in DC and X ! -!grad do k=i,j-1 -!grad do l=1,3 -!grad gvdwc(l,k)=gvdwc(l,k)+gg(l) -!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 + subroutine sc_grad_cat_pep + use calc_data + real(kind=8), dimension(3) :: dcosom1,dcosom2 + eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 & + +dCAVdOM1+ dGCLdOM1+ dPOLdOM1 + eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 & + +dCAVdOM2+ dGCLdOM2+ dPOLdOM2 + + eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 & + -2.0D0*alf12*eps3der+sigder*sigsq_om12& + +dCAVdOM12+ dGCLdOM12 +! diagnostics only +! eom1=0.0d0 +! eom2=0.0d0 +! eom12=evdwij*eps1_om12 +! end diagnostics + + do k=1,3 + dcosom1(k) = rij * (dc_norm(k,i) - om1 * erij(k)) + dcosom2(k) = rij * (dc_norm(k,nres+j) - om2 * erij(k)) + gg(k) = gg(k) + eom1 * dcosom1(k) + eom2 * dcosom2(k) + gvdwc_pepbase(k,i)= gvdwc_pepbase(k,i) +0.5*(- gg(k)) & + + (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i)))& + *dsci_inv*2.0 & + - (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0 + gvdwc_pepbase(k,i+1)= gvdwc_pepbase(k,i+1) +0.5*(- gg(k)) & + - (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i))) & + *dsci_inv*2.0 & + + (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0 + gradpepcat(k,j)=gradpepcat(k,j)+gg(k) + enddo + end subroutine sc_grad_cat_pep #ifdef CRYST_THETA !----------------------------------------------------------------------------- @@ -12633,7 +12651,7 @@ ! call intcartderiv ! call checkintcartgrad call zerograd - aincr=1.0D-4 + aincr=1.0D-5 write(iout,*) 'Calling CHECK_ECARTINT.' nf=0 icall=0 @@ -12860,7 +12878,7 @@ ! call intcartderiv ! call checkintcartgrad call zerograd - aincr=1.0D-7 + aincr=1.0D-6 write(iout,*) 'Calling CHECK_ECARTINT.',aincr nf=0 icall=0 @@ -22641,7 +22659,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, & @@ -22655,7 +22673,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 @@ -22667,7 +22685,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) @@ -22676,6 +22694,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) @@ -22752,7 +22772,7 @@ end subroutine ecatcat !--------------------------------------------------------------------------- ! new for K+ - subroutine ecats_prot_amber(ecations_prot_amber) + subroutine ecats_prot_amber(evdw) ! subroutine ecat_prot2(ecation_prot) use calc_data use comm_momo @@ -22782,7 +22802,7 @@ integer troll,jj,istate real (kind=8) :: dcosom1(3),dcosom2(3) - ecations_prot_amber=0.0D0 + evdw=0.0D0 if (nres_molec(5).eq.0) return eps_out=80.0d0 ! sss_ele_cut=1.0d0 @@ -22791,13 +22811,15 @@ do i=1,4 itmp=itmp+nres_molec(i) enddo + go to 17 ! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization do i=ibond_start,ibond_end ! print *,"I am in EVDW",i itypi=iabs(itype(i,1)) + ! if (i.ne.47) cycle - if (itypi.eq.ntyp1) cycle + if ((itypi.eq.ntyp1).or.(itypi.eq.10)) cycle itypi1=iabs(itype(i+1,1)) xi=c(1,nres+i) yi=c(2,nres+i) @@ -22815,11 +22837,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) @@ -22874,11 +22896,17 @@ ! alf2 = 0.0d0 ! alf12 = 0.0d0 ! a12sq = rborn(itypi,itypj) * rborn(itypj,itypi) - chi1 = chicat(itypi,itypj) - chis1 = chiscat(itypi,itypj) - chip1 = chippcat(itypi,itypj) + chi1 = chi1cat(itypi,itypj) + chis1 = chis1cat(itypi,itypj) + chip1 = chipp1cat(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 @@ -22972,11 +23000,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 @@ -23038,23 +23066,23 @@ END DO erdxi = scalar( ertail(1), dC_norm(1,i+nres) ) erdxj = scalar( ertail(1), dC_norm(1,j) ) - facd1 = dtail(1,itypi,itypj) * vbld_inv(i+nres) - facd2 = dtail(2,itypi,itypj) * vbld_inv(j+nres) + facd1 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres) + 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) + isel = iabs(Qi) + 1 ! ion is always charged so iabs(Qj) IF (isel.eq.0) THEN !c! No charges - do nothing eheadtail = 0.0d0 @@ -23090,7 +23118,7 @@ ELSE IF ((isel.eq.2.and. & iabs(Qi).eq.1).and. & - nstate(itypi,itypj).eq.1) THEN + nstatecat(itypi,itypj).eq.1) THEN !c! Same charge-charge interaction ( +/+ or -/- ) if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then @@ -23121,7 +23149,7 @@ ! ! CALL energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad) END IF ! this endif ends the "catch the gly-gly" at the beggining of Fcav - evdw = evdw + Fcav + eheadtail + evdw = evdw + Fcav + eheadtail IF (energy_dec) write (iout,'(2(1x,a3,i3),3f6.2,10f16.7)') & restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& @@ -23135,101 +23163,59 @@ !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. ! print *,"EVDW KURW",evdw,nres - - return - end subroutine ecats_prot_amber - -!--------------------------------------------------------------------------- -! old for Ca2+ - subroutine ecat_prot(ecation_prot) -! use calc_data -! use comm_momo - integer i,j,k,subchap,itmp,inum - real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,& - r7,r4,ecationcation - real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, & - dist_init,dist_temp,ecation_prot,rcal,rocal, & - Evan1,Evan2,EC,cm1mag,DASGL,delta,r0p,Epepcat, & - catl,cml,calpl, Etotal_p, Etotal_m,rtab,wdip,wmodquad,wquad1, & - wquad2,wvan1,E1,E2,wconst,wvan2,rcpm,dcmag,sin2thet,sinthet, & - costhet,v1m,v2m,wh2o,wc,rsecp,Ir,Irsecp,Irthrp,Irfourp,Irfiftp,& - Irsistp,Irseven,Irtwelv,Irthir,dE1dr,dE2dr,dEdcos,wquad2p,opt, & - rs,rthrp,rfourp,rsixp,reight,Irsixp,Ireight,Irtw,Irfourt, & - opt1,opt2,opt3,opt4,opt5,opt6,opt7,opt8,opt9,opt10,opt11,opt12,& - opt13,opt14,opt15,opt16,opt17,opt18,opt19, & - Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip,& - ndiv,ndivi - real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,& - gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, & - dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, & - tab1,tab2,tab3,diff,cm1,sc,p,tcat,talp,cm,drcp,drcp_norm,vcat, & - v1,v2,v3,myd_norm,dx,vcm,valpha,drdpep,dcosdpep,dcosddci,dEdpep,& - dEcCat,dEdipCm,dEdipCalp,dEquad1Cat,dEquad1Cm,dEquad1Calp, & - dEquad2Cat,dEquad2Cm,dEquad2Calpd,Evan1Cat,dEvan1Calp,dEvan2Cat,& - dEvan2Calp,dEtotalCat,dscvec,dEcCm,dEcCalp,dEdipCat,dEquad2Calp,& - dEvan1Cat - real(kind=8),dimension(6) :: vcatprm - ecation_prot=0.0d0 -! first lets calculate interaction with peptide groups - if (nres_molec(5).eq.0) return - itmp=0 - do i=1,4 - itmp=itmp+nres_molec(i) - enddo -! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization +!!! return + 17 continue do i=ibond_start,ibond_end -! cycle - if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle ! leave dummy atoms - xi=0.5d0*(c(1,i)+c(1,i+1)) - yi=0.5d0*(c(2,i)+c(2,i+1)) - zi=0.5d0*(c(3,i)+c(3,i+1)) - xi=mod(xi,boxxsize) + +! print *,"I am in EVDW",i + itypi=10 ! the peptide group parameters are for glicine + +! if (i.ne.47) cycle + if ((itype(i,1).eq.ntyp1).or.itype(i+1,1).eq.ntyp1) cycle + itypi1=iabs(itype(i+1,1)) + xi=(c(1,i)+c(1,i+1))/2.0 + yi=(c(2,i)+c(2,i+1))/2.0 + zi=(c(3,i)+c(3,i+1))/2.0 + xi=dmod(xi,boxxsize) if (xi.lt.0) xi=xi+boxxsize - yi=mod(yi,boxysize) + yi=dmod(yi,boxysize) if (yi.lt.0) yi=yi+boxysize - zi=mod(zi,boxzsize) + zi=dmod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize - + dxi=dc_norm(1,i) + dyi=dc_norm(2,i) + dzi=dc_norm(3,i) + dsci_inv=vbld_inv(i+1)/2.0 do j=itmp+1,itmp+nres_molec(5) -! print *,"WTF",itmp,j,i -! all parameters were for Ca2+ to approximate single charge divide by two - ndiv=1.0 - if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0 - wconst=78*ndiv - wdip =1.092777950857032D2 - wdip=wdip/wconst - wmodquad=-2.174122713004870D4 - wmodquad=wmodquad/wconst - wquad1 = 3.901232068562804D1 - wquad1=wquad1/wconst - wquad2 = 3 - wquad2=wquad2/wconst - wvan1 = 0.1 - wvan2 = 6 -! itmp=0 +! Calculate SC interaction energy. + itypj=iabs(itype(j,5)) + if ((itypj.eq.ntyp1)) cycle + CALL elgrad_init_cat_pep(eheadtail,Egb,Ecl,Elj,Equad,Epol) + + dscj_inv=0.0 xj=c(1,j) yj=c(2,j) zj=c(3,j) - xj=dmod(xj,boxxsize) - if (xj.lt.0) xj=xj+boxxsize - yj=dmod(yj,boxysize) - if (yj.lt.0) yj=yj+boxysize - zj=dmod(zj,boxzsize) - if (zj.lt.0) zj=zj+boxzsize - dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 - xj_safe=xj - yj_safe=yj - zj_safe=zj - subchap=0 - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 + xj=dmod(xj,boxxsize) + if (xj.lt.0) xj=xj+boxxsize + yj=dmod(yj,boxysize) + if (yj.lt.0) yj=yj+boxysize + zj=dmod(zj,boxzsize) + if (zj.lt.0) zj=zj+boxzsize + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 + + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 xj=xj_safe+xshift*boxxsize yj=yj_safe+yshift*boxysize zj=zj_safe+zshift*boxzsize @@ -23241,80 +23227,431 @@ zj_temp=zj subchap=1 endif - enddo - enddo - enddo - if (subchap.eq.1) then + enddo + enddo + enddo + if (subchap.eq.1) then xj=xj_temp-xi yj=yj_temp-yi zj=zj_temp-zi - else + else xj=xj_safe-xi yj=yj_safe-yi zj=zj_safe-zi - endif -! enddo -! enddo - rcpm = sqrt(xj**2+yj**2+zj**2) - drcp_norm(1)=xj/rcpm - drcp_norm(2)=yj/rcpm - drcp_norm(3)=zj/rcpm - dcmag=0.0 - do k=1,3 - dcmag=dcmag+dc(k,i)**2 - enddo - dcmag=dsqrt(dcmag) - do k=1,3 - myd_norm(k)=dc(k,i)/dcmag - enddo - costhet=drcp_norm(1)*myd_norm(1)+drcp_norm(2)*myd_norm(2)+& - drcp_norm(3)*myd_norm(3) - rsecp = rcpm**2 - Ir = 1.0d0/rcpm - Irsecp = 1.0d0/rsecp - Irthrp = Irsecp/rcpm - Irfourp = Irthrp/rcpm - Irfiftp = Irfourp/rcpm - Irsistp=Irfiftp/rcpm - Irseven=Irsistp/rcpm - Irtwelv=Irsistp*Irsistp - Irthir=Irtwelv/rcpm - sin2thet = (1-costhet*costhet) - sinthet=sqrt(sin2thet) - E1 = wdip*Irsecp*costhet+(wmodquad*Irfourp+wquad1*Irthrp)& - *sin2thet - E2 = -wquad1*Irthrp*wquad2+wvan1*(wvan2**12*Irtwelv-& - 2*wvan2**6*Irsistp) - ecation_prot = ecation_prot+E1+E2 -! print *,"ecatprot",i,j,ecation_prot,rcpm - dE1dr = -2*costhet*wdip*Irthrp-& - (4*wmodquad*Irfiftp+3*wquad1*Irfourp)*sin2thet - dE2dr = 3*wquad1*wquad2*Irfourp- & - 12*wvan1*wvan2**6*(wvan2**6*Irthir-Irseven) - dEdcos = wdip*Irsecp-2*(wmodquad*Irfourp+wquad1*Irthrp)*costhet - do k=1,3 - drdpep(k) = -drcp_norm(k) - dcosdpep(k) = Ir*(costhet*drcp_norm(k)-myd_norm(k)) - dcosddci(k) = drcp_norm(k)/dcmag-costhet*myd_norm(k)/dcmag - dEdpep(k) = (dE1dr+dE2dr)*drdpep(k)+dEdcos*dcosdpep(k) - dEddci(k) = dEdcos*dcosddci(k) - enddo - do k=1,3 - gradpepcat(k,i)=gradpepcat(k,i)+0.5D0*dEdpep(k)-dEddci(k) - gradpepcat(k,i+1)=gradpepcat(k,i+1)+0.5D0*dEdpep(k)+dEddci(k) - gradpepcat(k,j)=gradpepcat(k,j)-dEdpep(k) - enddo - enddo ! j - enddo ! i -!------------------------------------------sidechains -! do i=1,nres_molec(1) - do i=ibond_start,ibond_end - if ((itype(i,1).eq.ntyp1)) cycle ! leave dummy atoms -! cycle -! print *,i,ecation_prot - xi=(c(1,i+nres)) - yi=(c(2,i+nres)) - zi=(c(3,i+nres)) + endif + + dxj = 0.0d0! dc_norm( 1, nres+j ) + dyj = 0.0d0!dc_norm( 2, nres+j ) + dzj = 0.0d0! dc_norm( 3, nres+j ) + + itypi = 10 + itypj = itype(j,5) +! Parameters from fitting the analitical expressions to the PMF obtained by umbrella +! sampling performed with amber package +! alf1 = 0.0d0 +! alf2 = 0.0d0 +! alf12 = 0.0d0 +! a12sq = rborn(itypi,itypj) * rborn(itypj,itypi) + chi1 = chi1cat(itypi,itypj) + chis1 = chis1cat(itypi,itypj) + chip1 = chipp1cat(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 + sig1 = sigmap1cat(itypi,itypj) +! sig2 = sigmap2(itypi,itypj) +! alpha factors from Fcav/Gcav + b1cav = alphasurcat(1,itypi,itypj) + b2cav = alphasurcat(2,itypi,itypj) + b3cav = alphasurcat(3,itypi,itypj) + b4cav = alphasurcat(4,itypi,itypj) + +! used to determine whether we want to do quadrupole calculations + eps_in = epsintabcat(itypi,itypj) + if (eps_in.eq.0.0) eps_in=1.0 + + eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out)) +! Rtail = 0.0d0 + + DO k = 1, 3 + ctail(k,1)=(c(k,i)+c(k,i+1))/2.0 + ctail(k,2)=c(k,j) + END DO +!c! tail distances will be themselves usefull elswhere +!c1 (in Gcav, for example) + Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 ) + Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 ) + Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 ) + Rtail = dsqrt( & + (Rtail_distance(1)*Rtail_distance(1)) & + + (Rtail_distance(2)*Rtail_distance(2)) & + + (Rtail_distance(3)*Rtail_distance(3))) +! tail location and distance calculations +! dhead1 + d1 = dheadcat(1, 1, itypi, itypj) +! d2 = dhead(2, 1, itypi, itypj) + DO k = 1,3 +! location of polar head is computed by taking hydrophobic centre +! and moving by a d1 * dc_norm vector +! see unres publications for very informative images + chead(k,1) = (c(k, i)+c(k,i+1))/2.0 + d1 * dc_norm(k, i) + chead(k,2) = c(k, j) +! distance +! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres)) +! Rsc(k) = Rsc_distance(k) * Rsc_distance(k) + Rhead_distance(k) = chead(k,2) - chead(k,1) + END DO +! pitagoras (root of sum of squares) + Rhead = dsqrt( & + (Rhead_distance(1)*Rhead_distance(1)) & + + (Rhead_distance(2)*Rhead_distance(2)) & + + (Rhead_distance(3)*Rhead_distance(3))) +!------------------------------------------------------------------- +! zero everything that should be zero'ed + evdwij = 0.0d0 + ECL = 0.0d0 + Elj = 0.0d0 + Equad = 0.0d0 + Epol = 0.0d0 + Fcav=0.0d0 + eheadtail = 0.0d0 + dGCLdOM1 = 0.0d0 + dGCLdOM2 = 0.0d0 + dGCLdOM12 = 0.0d0 + dPOLdOM1 = 0.0d0 + dPOLdOM2 = 0.0d0 + Fcav = 0.0d0 + dFdR = 0.0d0 + dCAVdOM1 = 0.0d0 + dCAVdOM2 = 0.0d0 + dCAVdOM12 = 0.0d0 + dscj_inv = vbld_inv(j+nres) +! print *,i,j,dscj_inv,dsci_inv +! rij holds 1/(distance of Calpha atoms) + rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj) + rij = dsqrt(rrij) + CALL sc_angular +! this should be in elgrad_init but om's are calculated by sc_angular +! which in turn is used by older potentials +! om = omega, sqom = om^2 + sqom1 = om1 * om1 + sqom2 = om2 * om2 + sqom12 = om12 * om12 + +! now we calculate EGB - Gey-Berne +! It will be summed up in evdwij and saved in evdw + sigsq = 1.0D0 / sigsq + sig = sig0ij * dsqrt(sigsq) +! rij_shift = 1.0D0 / rij - sig + sig0ij + rij_shift = Rtail - sig + sig0ij + IF (rij_shift.le.0.0D0) THEN + evdw = 1.0D20 + RETURN + END IF + sigder = -sig * sigsq + rij_shift = 1.0D0 / rij_shift + fac = rij_shift**expon + c1 = fac * fac * aa_aq_cat(itypi,itypj) +! print *,"ADAM",aa_aq(itypi,itypj) + +! c1 = 0.0d0 + c2 = fac * bb_aq_cat(itypi,itypj) +! c2 = 0.0d0 + evdwij = eps1 * eps2rt * eps3rt * ( c1 + c2 ) + eps2der = eps3rt * evdwij + eps3der = eps2rt * evdwij +! evdwij = 4.0d0 * eps2rt * eps3rt * evdwij + evdwij = eps2rt * eps3rt * evdwij +!#ifdef TSCSC +! IF (bb_aq(itypi,itypj).gt.0) THEN +! evdw_p = evdw_p + evdwij +! ELSE +! evdw_m = evdw_m + evdwij +! END IF +!#else + evdw = evdw & + + evdwij +!#endif + c1 = c1 * eps1 * eps2rt**2 * eps3rt**2 + fac = -expon * (c1 + evdwij) * rij_shift + sigder = fac * sigder +! Calculate distance derivative + gg(1) = fac + gg(2) = fac + gg(3) = fac + + fac = chis1 * sqom1 + chis2 * sqom2 & + - 2.0d0 * chis12 * om1 * om2 * om12 + + pom = 1.0d0 - chis1 * chis2 * sqom12 + print *,"TUT2",fac,chis1,sqom1,pom + Lambf = (1.0d0 - (fac / pom)) + Lambf = dsqrt(Lambf) + sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0) + Chif = Rtail * sparrow + ChiLambf = Chif * Lambf + eagle = dsqrt(ChiLambf) + bat = ChiLambf ** 11.0d0 + top = b1cav * ( eagle + b2cav * ChiLambf - b3cav ) + bot = 1.0d0 + b4cav * (ChiLambf ** 12.0d0) + botsq = bot * bot + Fcav = top / bot + + dtop = b1cav * ((Lambf / (2.0d0 * eagle)) + (b2cav * Lambf)) + dbot = 12.0d0 * b4cav * bat * Lambf + dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow + + dtop = b1cav * ((Chif / (2.0d0 * eagle)) + (b2cav * Chif)) + dbot = 12.0d0 * b4cav * bat * Chif + eagle = Lambf * pom + dFdOM1 = -(chis1 * om1 - chis12 * om2 * om12) / (eagle) + dFdOM2 = -(chis2 * om2 - chis12 * om1 * om12) / (eagle) + dFdOM12 = chis12 * (chis1 * om1 * om12 - om2) & + * (chis2 * om2 * om12 - om1) / (eagle * pom) + + dFdL = ((dtop * bot - top * dbot) / botsq) + dCAVdOM1 = dFdL * ( dFdOM1 ) + dCAVdOM2 = dFdL * ( dFdOM2 ) + dCAVdOM12 = dFdL * ( dFdOM12 ) + + DO k= 1, 3 + ertail(k) = Rtail_distance(k)/Rtail + END DO + erdxi = scalar( ertail(1), dC_norm(1,i) ) + erdxj = scalar( ertail(1), dC_norm(1,j) ) + facd1 = dtailcat(1,itypi,itypj) * vbld_inv(i) + facd2 = dtailcat(2,itypi,itypj) * vbld_inv(j+nres) + DO k = 1, 3 + pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(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) + gradpepcat(k,i) = gradpepcat(k,i) & + - (( dFdR + gg(k) ) * ertail(k))/2.0d0 + gradpepcat(k,i+1) = gradpepcat(k,i+1) & + - (( dFdR + gg(k) ) * ertail(k))/2.0d0 + + 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 = 3 +!c! Dipole-charge interactions + if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then + Qi=Qi*2 + Qij=Qij*2 + endif + if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then + Qj=Qj*2 + Qij=Qij*2 + endif + CALL edq_cat_pep(ecl, elj, epol) + eheadtail = ECL + elj + epol +! print *,"i,",i,eheadtail + eheadtail = 0.0d0 + + evdw = evdw + Fcav + eheadtail + + IF (energy_dec) write (iout,'(2(1x,a3,i3),3f6.2,10f16.7)') & + restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& + 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& + Equad,evdwij+Fcav+eheadtail,evdw +! evdw = evdw + Fcav + eheadtail + +! iF (nstate(itypi,itypj).eq.1) THEN + CALL sc_grad_cat_pep +! END IF +!c!------------------------------------------------------------------- +!c! NAPISY KONCOWE + END DO ! j + END DO ! i +!c write (iout,*) "Number of loop steps in EGB:",ind +!c energy_dec=.false. +! print *,"EVDW KURW",evdw,nres + + + return + end subroutine ecats_prot_amber + +!--------------------------------------------------------------------------- +! old for Ca2+ + subroutine ecat_prot(ecation_prot) +! use calc_data +! use comm_momo + integer i,j,k,subchap,itmp,inum + real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,& + r7,r4,ecationcation + real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, & + dist_init,dist_temp,ecation_prot,rcal,rocal, & + Evan1,Evan2,EC,cm1mag,DASGL,delta,r0p,Epepcat, & + catl,cml,calpl, Etotal_p, Etotal_m,rtab,wdip,wmodquad,wquad1, & + wquad2,wvan1,E1,E2,wconst,wvan2,rcpm,dcmag,sin2thet,sinthet, & + costhet,v1m,v2m,wh2o,wc,rsecp,Ir,Irsecp,Irthrp,Irfourp,Irfiftp,& + Irsistp,Irseven,Irtwelv,Irthir,dE1dr,dE2dr,dEdcos,wquad2p,opt, & + rs,rthrp,rfourp,rsixp,reight,Irsixp,Ireight,Irtw,Irfourt, & + opt1,opt2,opt3,opt4,opt5,opt6,opt7,opt8,opt9,opt10,opt11,opt12,& + opt13,opt14,opt15,opt16,opt17,opt18,opt19, & + Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip,& + ndiv,ndivi + real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,& + gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, & + dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, & + tab1,tab2,tab3,diff,cm1,sc,p,tcat,talp,cm,drcp,drcp_norm,vcat, & + v1,v2,v3,myd_norm,dx,vcm,valpha,drdpep,dcosdpep,dcosddci,dEdpep,& + dEcCat,dEdipCm,dEdipCalp,dEquad1Cat,dEquad1Cm,dEquad1Calp, & + dEquad2Cat,dEquad2Cm,dEquad2Calpd,Evan1Cat,dEvan1Calp,dEvan2Cat,& + dEvan2Calp,dEtotalCat,dscvec,dEcCm,dEcCalp,dEdipCat,dEquad2Calp,& + dEvan1Cat + real(kind=8),dimension(6) :: vcatprm + ecation_prot=0.0d0 +! first lets calculate interaction with peptide groups + if (nres_molec(5).eq.0) return + itmp=0 + do i=1,4 + itmp=itmp+nres_molec(i) + enddo +! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization + do i=ibond_start,ibond_end +! cycle + if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle ! leave dummy atoms + xi=0.5d0*(c(1,i)+c(1,i+1)) + yi=0.5d0*(c(2,i)+c(2,i+1)) + zi=0.5d0*(c(3,i)+c(3,i+1)) + xi=mod(xi,boxxsize) + if (xi.lt.0) xi=xi+boxxsize + yi=mod(yi,boxysize) + if (yi.lt.0) yi=yi+boxysize + zi=mod(zi,boxzsize) + if (zi.lt.0) zi=zi+boxzsize + + do j=itmp+1,itmp+nres_molec(5) +! print *,"WTF",itmp,j,i +! all parameters were for Ca2+ to approximate single charge divide by two + ndiv=1.0 + if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0 + wconst=78*ndiv + wdip =1.092777950857032D2 + wdip=wdip/wconst + wmodquad=-2.174122713004870D4 + wmodquad=wmodquad/wconst + wquad1 = 3.901232068562804D1 + wquad1=wquad1/wconst + wquad2 = 3 + wquad2=wquad2/wconst + wvan1 = 0.1 + wvan2 = 6 +! itmp=0 + + xj=c(1,j) + yj=c(2,j) + zj=c(3,j) + xj=dmod(xj,boxxsize) + if (xj.lt.0) xj=xj+boxxsize + yj=dmod(yj,boxysize) + if (yj.lt.0) yj=yj+boxysize + zj=dmod(zj,boxzsize) + if (zj.lt.0) zj=zj+boxzsize + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif +! enddo +! enddo + rcpm = sqrt(xj**2+yj**2+zj**2) + drcp_norm(1)=xj/rcpm + drcp_norm(2)=yj/rcpm + drcp_norm(3)=zj/rcpm + dcmag=0.0 + do k=1,3 + dcmag=dcmag+dc(k,i)**2 + enddo + dcmag=dsqrt(dcmag) + do k=1,3 + myd_norm(k)=dc(k,i)/dcmag + enddo + costhet=drcp_norm(1)*myd_norm(1)+drcp_norm(2)*myd_norm(2)+& + drcp_norm(3)*myd_norm(3) + rsecp = rcpm**2 + Ir = 1.0d0/rcpm + Irsecp = 1.0d0/rsecp + Irthrp = Irsecp/rcpm + Irfourp = Irthrp/rcpm + Irfiftp = Irfourp/rcpm + Irsistp=Irfiftp/rcpm + Irseven=Irsistp/rcpm + Irtwelv=Irsistp*Irsistp + Irthir=Irtwelv/rcpm + sin2thet = (1-costhet*costhet) + sinthet=sqrt(sin2thet) + E1 = wdip*Irsecp*costhet+(wmodquad*Irfourp+wquad1*Irthrp)& + *sin2thet + E2 = -wquad1*Irthrp*wquad2+wvan1*(wvan2**12*Irtwelv-& + 2*wvan2**6*Irsistp) + ecation_prot = ecation_prot+E1+E2 +! print *,"ecatprot",i,j,ecation_prot,rcpm + dE1dr = -2*costhet*wdip*Irthrp-& + (4*wmodquad*Irfiftp+3*wquad1*Irfourp)*sin2thet + dE2dr = 3*wquad1*wquad2*Irfourp- & + 12*wvan1*wvan2**6*(wvan2**6*Irthir-Irseven) + dEdcos = wdip*Irsecp-2*(wmodquad*Irfourp+wquad1*Irthrp)*costhet + do k=1,3 + drdpep(k) = -drcp_norm(k) + dcosdpep(k) = Ir*(costhet*drcp_norm(k)-myd_norm(k)) + dcosddci(k) = drcp_norm(k)/dcmag-costhet*myd_norm(k)/dcmag + dEdpep(k) = (dE1dr+dE2dr)*drdpep(k)+dEdcos*dcosdpep(k) + dEddci(k) = dEdcos*dcosddci(k) + enddo + do k=1,3 + gradpepcat(k,i)=gradpepcat(k,i)+0.5D0*dEdpep(k)-dEddci(k) + gradpepcat(k,i+1)=gradpepcat(k,i+1)+0.5D0*dEdpep(k)+dEddci(k) + gradpepcat(k,j)=gradpepcat(k,j)-dEdpep(k) + enddo + enddo ! j + enddo ! i +!------------------------------------------sidechains +! do i=1,nres_molec(1) + do i=ibond_start,ibond_end + if ((itype(i,1).eq.ntyp1)) cycle ! leave dummy atoms +! cycle +! print *,i,ecation_prot + xi=(c(1,i+nres)) + yi=(c(2,i+nres)) + zi=(c(3,i+nres)) xi=mod(xi,boxxsize) if (xi.lt.0) xi=xi+boxxsize yi=mod(yi,boxysize) @@ -26347,7 +26684,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& @@ -26357,13 +26694,13 @@ - dGLJdR * pom pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j)) - gvdwx(k,j) = gvdwx(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 +! 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)& @@ -26371,7 +26708,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) & @@ -26903,15 +27240,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) & - + dPOLdR2 * condor +! 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 @@ -27059,8 +27396,8 @@ !c!------------------------------------------------------------------- !c! ecl - sparrow = w1 * Qi * om1 - hawk = w2 * Qi * Qi * (1.0d0 - sqom2) + sparrow = w1 * Qj * om1 + hawk = w2 * Qj * Qj * (1.0d0 - sqom2) ECL = sparrow / Rhead**2.0d0 & - hawk / Rhead**4.0d0 !c!------------------------------------------------------------------- @@ -27069,9 +27406,9 @@ dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 & + 4.0d0 * hawk / Rhead**5.0d0 !c! dF/dom1 - dGCLdOM1 = (w1 * Qi) / (Rhead**2.0d0) + dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0) !c! dF/dom2 - dGCLdOM2 = (2.0d0 * w2 * Qi * Qi * om2) / (Rhead ** 4.0d0) + dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0) !c-------------------------------------------------------------------- !c Polarization energy !c Epol @@ -27175,8 +27512,8 @@ !c!------------------------------------------------------------------- !c! ecl - sparrow = w1 * Qi * om1 - hawk = w2 * Qi * Qi * (1.0d0 - sqom2) + sparrow = w1 * Qj * om1 + hawk = w2 * Qj * Qj * (1.0d0 - sqom2) ECL = sparrow / Rhead**2.0d0 & - hawk / Rhead**4.0d0 !c!------------------------------------------------------------------- @@ -27185,9 +27522,9 @@ dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 & + 4.0d0 * hawk / Rhead**5.0d0 !c! dF/dom1 - dGCLdOM1 = (w1 * Qi) / (Rhead**2.0d0) + dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0) !c! dF/dom2 - dGCLdOM2 = (2.0d0 * w2 * Qi * Qi * om2) / (Rhead ** 4.0d0) + dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0) !c-------------------------------------------------------------------- !c-------------------------------------------------------------------- !c Polarization energy @@ -27232,31 +27569,31 @@ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) ) facd1 = d1 * vbld_inv(i+nres) facd2 = d2 * vbld_inv(j) - facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres) + facd3 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres) DO k = 1, 3 condor = (erhead_tail(k,2) & + 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) & - + dGCLdR * pom & - + dPOLdR2 * condor & - + dGLJdR * pom +! 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) @@ -27265,6 +27602,129 @@ RETURN END SUBROUTINE edq_cat + SUBROUTINE edq_cat_pep(Ecl,Elj,Epol) + use comm_momo + use calc_data + + double precision facd3, adler,ecl,elj,epol + alphapol2 = alphapolcat(itypj,itypi) + w1 = wqdipcat(1,itypi,itypj) + w2 = wqdipcat(2,itypi,itypj) + pis = sig0headcat(itypi,itypj) + eps_head = epsheadcat(itypi,itypj) +!c!------------------------------------------------------------------- +!c! R2 - distance between head of jth side chain and tail of ith sidechain + R2 = 0.0d0 + DO k = 1, 3 +!c! Calculate head-to-tail distances + R2=R2+(chead(k,2)-ctail(k,1))**2 + END DO +!c! Pitagoras + R2 = dsqrt(R2) + +!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj) +!c! & +dhead(1,1,itypi,itypj))**2)) +!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj) +!c! & +dhead(2,1,itypi,itypj))**2)) + + +!c!------------------------------------------------------------------- +!c! ecl + sparrow = w1 * Qj * om1 + hawk = w2 * Qj * Qj * (1.0d0 - sqom2) +! print *,"CO?!.", w1,w2,Qj,om1 + ECL = sparrow / Rhead**2.0d0 & + - hawk / Rhead**4.0d0 +!c!------------------------------------------------------------------- +!c! derivative of ecl is Gcl +!c! dF/dr part + dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 & + + 4.0d0 * hawk / Rhead**5.0d0 +!c! dF/dom1 + dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0) +!c! dF/dom2 + dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0) +!c-------------------------------------------------------------------- +!c-------------------------------------------------------------------- +!c Polarization energy +!c Epol + MomoFac2 = (1.0d0 - chi2 * sqom1) + RR2 = R2 * R2 / MomoFac2 + ee2 = exp(-(RR2 / (4.0d0 * a12sq))) + fgb2 = sqrt(RR2 + a12sq * ee2) + epol = 332.0d0 * eps_inout_fac * ((alphapol2/fgb2) ** 4.0d0 ) + dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0) & + / (fgb2 ** 5.0d0) + dFGBdR2 = ( (R2 / MomoFac2) & + * ( 2.0d0 - (0.5d0 * ee2) ) ) & + / (2.0d0 * fgb2) + dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) & + * (2.0d0 - 0.5d0 * ee2) ) & + / (2.0d0 * fgb2) + dPOLdR2 = dPOLdFGB2 * dFGBdR2 +!c! dPOLdR2 = 0.0d0 + dPOLdOM1 = dPOLdFGB2 * dFGBdOM1 +!c! dPOLdOM1 = 0.0d0 + dPOLdOM2 = 0.0d0 +!c!------------------------------------------------------------------- +!c! Elj + pom = (pis / Rhead)**6.0d0 + Elj = 4.0d0 * eps_head * pom * (pom-1.0d0) +!c! derivative of Elj is Glj + dGLJdR = 4.0d0 * eps_head & + * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) & + + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0))) +!c!------------------------------------------------------------------- + +!c! Return the results +!c! (see comments in Eqq) + DO k = 1, 3 + erhead(k) = Rhead_distance(k)/Rhead + erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2) + END DO + erdxi = scalar( erhead(1), dC_norm(1,i+nres) ) + erdxj = scalar( erhead(1), dC_norm(1,j) ) + eagle = scalar( erhead_tail(1,2), dC_norm(1,j) ) + adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) ) + facd1 = d1 * vbld_inv(i+1)/2.0 + facd2 = d2 * vbld_inv(j) + facd3 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres) + DO k = 1, 3 + condor = (erhead_tail(k,2) & + + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j))) + + pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres)) +! 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)) +! gradpepcatx(k,j) = gradpepcatx(k,j) & +! + dGCLdR * pom & +! + dPOLdR2 * condor & +! + dGLJdR * pom + + + gradpepcat(k,i) = gradpepcat(k,i) +0.5d0*( & + - dGCLdR * erhead(k) & + - dPOLdR2 * erhead_tail(k,2) & + - dGLJdR * erhead(k)) + gradpepcat(k,i+1) = gradpepcat(k,i+1) +0.5d0*( & + - dGCLdR * erhead(k) & + - dPOLdR2 * erhead_tail(k,2) & + - dGLJdR * erhead(k)) + + + gradpepcat(k,j) = gradpepcat(k,j) & + + dGCLdR * erhead(k) & + + dPOLdR2 * erhead_tail(k,2) & + + dGLJdR * erhead(k) + + END DO + RETURN + END SUBROUTINE edq_cat_pep SUBROUTINE edd(ECL) ! IMPLICIT NONE @@ -27472,7 +27932,7 @@ real(kind=8) :: eheadtail,Egb,Ecl,Elj,Equad,Epol,Rb eps_out=80.0d0 itypi = itype(i,1) - itypj = itype(j,1) + itypj = itype(j,5) !c! 1/(Gas Constant * Thermostate temperature) = BetaT !c! ENABLE THIS LINE WHEN USING CHECKGRAD!!! !c! t_bath = 300 @@ -27481,31 +27941,17 @@ BetaT = 1.0d0 / (298.0d0 * Rb) !c! Gay-berne var's sig0ij = sigmacat( itypi,itypj ) - chi1 = chicat( itypi, itypj ) -! chi2 = chi( itypj, itypi ) + chi1 = chi1cat( itypi, itypj ) chi2 = 0.0d0 -! chi12 = chi1 * chi2 chi12 = 0.0d0 - chip1 = chippcat( itypi, itypj ) -! chip2 = chipp( itypj, itypi ) + chip1 = chipp1cat( itypi, itypj ) chip2 = 0.0d0 -! chip12 = chip1 * chip2 chip12 = 0.0d0 -! chi1=0.0 -! chi2=0.0 -! chi12=0.0 -! chip1=0.0 -! chip2=0.0 -! chip12=0.0 !c! not used by momo potential, but needed by sc_angular which is shared !c! by all energy_potential subroutines alf1 = 0.0d0 alf2 = 0.0d0 alf12 = 0.0d0 -!c! location, location, location -! xj = c( 1, nres+j ) - xi -! yj = c( 2, nres+j ) - yi -! zj = c( 3, nres+j ) - zi dxj = dc_norm( 1, nres+j ) dyj = dc_norm( 2, nres+j ) dzj = dc_norm( 3, nres+j ) @@ -27513,23 +27959,20 @@ d1 = dheadcat(1, 1, itypi, itypj) d2 = dheadcat(2, 1, itypi, itypj) !c! ai*aj from Fgb - a12sq = rborncat(itypi,itypj) * rborncat(itypj,itypi) + a12sq = rborn1cat(itypi,itypj) * rborn2cat(itypi,itypj) !c! a12sq = a12sq * a12sq !c! charge of amino acid itypi is... - Qi = ichargecat(itypi) + Qi = icharge(itypi) Qj = ichargecat(itypj) Qij = Qi * Qj !c! chis1,2,12 - chis1 = chiscat(itypi,itypj) -! chis2 = chis(itypj,itypi) + chis1 = chis1cat(itypi,itypj) chis2 = 0.0d0 -! chis12 = chis1 * chis2 chis12 = 0.0d0 sig1 = sigmap1cat(itypi,itypj) sig2 = sigmap2cat(itypi,itypj) !c! alpha factors from Fcav/Gcav b1cav = alphasurcat(1,itypi,itypj) -! b1cav=0.0 b2cav = alphasurcat(2,itypi,itypj) b3cav = alphasurcat(3,itypi,itypj) b4cav = alphasurcat(4,itypi,itypj) @@ -27542,7 +27985,7 @@ Rtail = 0.0d0 DO k = 1, 3 ctail(k,1)=c(k,i+nres)-dtailcat(1,itypi,itypj)*dc_norm(k,nres+i) - ctail(k,2)=c(k,j+nres)-dtailcat(2,itypi,itypj)*dc_norm(k,nres+j) + ctail(k,2)=c(k,j)!-dtailcat(2,itypi,itypj)*dc_norm(k,nres+j) END DO !c! tail distances will be themselves usefull elswhere !c1 (in Gcav, for example) @@ -27565,7 +28008,7 @@ !c! and moving by a d1 * dc_norm vector !c! see unres publications for very informative images chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres) - chead(k,2) = c(k, j+nres) + d2 * dc_norm(k, j+nres) + chead(k,2) = c(k, j) !c! distance !c! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres)) !c! Rsc(k) = Rsc_distance(k) * Rsc_distance(k) @@ -27592,6 +28035,114 @@ RETURN END SUBROUTINE elgrad_init_cat + SUBROUTINE elgrad_init_cat_pep(eheadtail,Egb,Ecl,Elj,Equad,Epol) + use comm_momo + use calc_data + real(kind=8) :: eheadtail,Egb,Ecl,Elj,Equad,Epol,Rb + eps_out=80.0d0 + itypi = 10 + itypj = itype(j,5) +!c! 1/(Gas Constant * Thermostate temperature) = BetaT +!c! ENABLE THIS LINE WHEN USING CHECKGRAD!!! +!c! t_bath = 300 +!c! BetaT = 1.0d0 / (t_bath * Rb)i + Rb=0.001986d0 + BetaT = 1.0d0 / (298.0d0 * Rb) +!c! Gay-berne var's + sig0ij = sigmacat( itypi,itypj ) + chi1 = chi1cat( itypi, itypj ) + chi2 = 0.0d0 + chi12 = 0.0d0 + chip1 = chipp1cat( itypi, itypj ) + chip2 = 0.0d0 + chip12 = 0.0d0 +!c! not used by momo potential, but needed by sc_angular which is shared +!c! by all energy_potential subroutines + alf1 = 0.0d0 + alf2 = 0.0d0 + alf12 = 0.0d0 + dxj = 0.0d0 !dc_norm( 1, nres+j ) + dyj = 0.0d0 !dc_norm( 2, nres+j ) + dzj = 0.0d0 !dc_norm( 3, nres+j ) +!c! distance from center of chain(?) to polar/charged head + d1 = dheadcat(1, 1, itypi, itypj) + d2 = dheadcat(2, 1, itypi, itypj) +!c! ai*aj from Fgb + a12sq = rborn1cat(itypi,itypj) * rborn2cat(itypi,itypj) +!c! a12sq = a12sq * a12sq +!c! charge of amino acid itypi is... + Qi = 0 + Qj = ichargecat(itypj) +! Qij = Qi * Qj +!c! chis1,2,12 + chis1 = chis1cat(itypi,itypj) + chis2 = 0.0d0 + chis12 = 0.0d0 + sig1 = sigmap1cat(itypi,itypj) + sig2 = sigmap2cat(itypi,itypj) +!c! alpha factors from Fcav/Gcav + b1cav = alphasurcat(1,itypi,itypj) + b2cav = alphasurcat(2,itypi,itypj) + b3cav = alphasurcat(3,itypi,itypj) + b4cav = alphasurcat(4,itypi,itypj) + wqd = wquadcat(itypi, itypj) +!c! used by Fgb + eps_in = epsintabcat(itypi,itypj) + eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out)) +!c!------------------------------------------------------------------- +!c! tail location and distance calculations + Rtail = 0.0d0 + DO k = 1, 3 + ctail(k,1)=(c(k,i)+c(k,i+1))/2.0-dtailcat(1,itypi,itypj)*dc_norm(k,i) + ctail(k,2)=c(k,j)!-dtailcat(2,itypi,itypj)*dc_norm(k,nres+j) + END DO +!c! tail distances will be themselves usefull elswhere +!c1 (in Gcav, for example) + Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 ) + Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 ) + Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 ) + Rtail = dsqrt( & + (Rtail_distance(1)*Rtail_distance(1)) & + + (Rtail_distance(2)*Rtail_distance(2)) & + + (Rtail_distance(3)*Rtail_distance(3))) +!c!------------------------------------------------------------------- +!c! Calculate location and distance between polar heads +!c! distance between heads +!c! for each one of our three dimensional space... + d1 = dheadcat(1, 1, itypi, itypj) + d2 = dheadcat(2, 1, itypi, itypj) + + DO k = 1,3 +!c! location of polar head is computed by taking hydrophobic centre +!c! and moving by a d1 * dc_norm vector +!c! see unres publications for very informative images + chead(k,1) = (c(k, i)+c(k,i+1))/2.0 + d1 * dc_norm(k, i) + chead(k,2) = c(k, j) +!c! distance +!c! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres)) +!c! Rsc(k) = Rsc_distance(k) * Rsc_distance(k) + Rhead_distance(k) = chead(k,2) - chead(k,1) + END DO +!c! pitagoras (root of sum of squares) + Rhead = dsqrt( & + (Rhead_distance(1)*Rhead_distance(1)) & + + (Rhead_distance(2)*Rhead_distance(2)) & + + (Rhead_distance(3)*Rhead_distance(3))) +!c!------------------------------------------------------------------- +!c! zero everything that should be zero'ed + Egb = 0.0d0 + ECL = 0.0d0 + Elj = 0.0d0 + Equad = 0.0d0 + Epol = 0.0d0 + eheadtail = 0.0d0 + dGCLdOM1 = 0.0d0 + dGCLdOM2 = 0.0d0 + dGCLdOM12 = 0.0d0 + dPOLdOM1 = 0.0d0 + dPOLdOM2 = 0.0d0 + RETURN + END SUBROUTINE elgrad_init_cat_pep double precision function tschebyshev(m,n,x,y) implicit none