X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.f90;h=d77030768df48ac79067415e8de7abc3b3a68ee1;hb=4baa9e481f53d4c89b076f3aece756fc47282649;hp=964624f7487c440f280a766c0bbdfdcc6f1fee29;hpb=8ca97b16fe25b7053f258263899ba030572cc58f;p=unres4.git diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index 964624f..d770307 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -130,9 +130,16 @@ gvdwpsb1,gelpp,gvdwpsb,gelsbc,gelsbx,gvdwsbx,gvdwsbc,gsbloc,& gsblocx,gradcorr_nucl,gradxorr_nucl,gradcorr3_nucl,gradxorr3_nucl,& gvdwpp_nucl +!-----------------------------NUCLEIC-PROTEIN GRADIENT + real(kind=8),dimension(:,:),allocatable :: gvdwx_scbase,gvdwc_scbase,& + gvdwx_pepbase,gvdwc_pepbase,gvdwx_scpho,gvdwc_scpho,& + gvdwc_peppho !------------------------------IONS GRADIENT - real(kind=8),dimension(:,:),allocatable :: gradcatcat + real(kind=8),dimension(:,:),allocatable :: gradcatcat, & + gradpepcat,gradpepcatx ! real(kind=8),dimension(:,:),allocatable :: gloc,gloc_x !(maxvar,2) + + real(kind=8),dimension(:,:),allocatable :: gel_loc,gel_loc_long,& gcorr3_turn,gcorr4_turn,gcorr6_turn,gradb,gradbx !(3,maxres) real(kind=8),dimension(:),allocatable :: gel_loc_loc,& @@ -238,6 +245,11 @@ real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,& ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,& ecorr3_nucl +! energies for ions + real(kind=8) :: ecation_prot,ecationcation +! energies for protein nucleic acid interaction + real(kind=8) :: escbase,epepbase,escpho,epeppho + #ifdef MPI real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw ! shielding effect varibles for MPI @@ -290,6 +302,13 @@ weights_(36)=wtor_d_nucl weights_(37)=wcorr_nucl weights_(38)=wcorr3_nucl + weights_(41)=wcatcat + weights_(42)=wcatprot + weights_(46)=wscbase + weights_(47)=wscpho + weights_(48)=wpeppho +! wcatcat= weights(41) +! wcatprot=weights(42) ! FG Master broadcasts the WEIGHTS_ array call MPI_Bcast(weights_(1),n_ene,& @@ -330,7 +349,11 @@ wtor_d_nucl =weights(36) wcorr_nucl =weights(37) wcorr3_nucl =weights(38) - + wcatcat= weights(41) + wcatprot=weights(42) + wscbase=weights(46) + wscpho=weights(47) + wpeppho=weights(48) endif time_Bcast=time_Bcast+MPI_Wtime()-time00 time_Bcastw=time_Bcastw+MPI_Wtime()-time00 @@ -581,7 +604,9 @@ etube=0.0d0 endif !-------------------------------------------------------- +! write (iout,*) "NRES_MOLEC(2),",nres_molec(2) ! print *,"before",ees,evdw1,ecorr + if (nres_molec(2).gt.0) then call ebond_nucl(estr_nucl) call ebend_nucl(ebe_nucl) call etor_nucl(etors_nucl) @@ -590,7 +615,29 @@ call epsb(evdwpsb,eelpsb) call esb(esbloc) call multibody_hb_nucl(ecorr_nucl,ecorr3_nucl,n_corr,n_corr1) - + else + etors_nucl=0.0d0 + estr_nucl=0.0d0 + ebe_nucl=0.0d0 + evdwsb=0.0d0 + eelsb=0.0d0 + esbloc=0.0d0 + endif + if (nfgtasks.gt.1) then + if (fg_rank.eq.0) then + call ecatcat(ecationcation) + endif + else + call ecatcat(ecationcation) + endif + call ecat_prot(ecation_prot) + if (nres_molec(2).gt.0) then + call eprot_sc_base(escbase) + call epep_sc_base(epepbase) + call eprot_sc_phosphate(escpho) + call eprot_pep_phosphate(epeppho) + endif +! call ecatcat(ecationcation) ! print *,"after ebend", ebe_nucl #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 @@ -655,6 +702,12 @@ ! Here are the energies showed per procesor if the are more processors ! per molecule then we sum it up in sum_energy subroutine ! print *," Processor",myrank," calls SUM_ENERGY" + energia(41)=ecation_prot + energia(42)=ecationcation + energia(46)=escbase + energia(47)=epepbase + energia(48)=escpho + energia(49)=epeppho call sum_energy(energia,.true.) if (dyn_ss) call dyn_set_nss ! print *," Processor",myrank," left SUM_ENERGY" @@ -697,7 +750,8 @@ real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,& ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,& ecorr3_nucl - + real(kind=8) :: ecation_prot,ecationcation + real(kind=8) :: escbase,epepbase,escpho,epeppho integer :: i #ifdef MPI integer :: ierr @@ -774,6 +828,14 @@ etors_d_nucl=energia(36) ecorr_nucl=energia(37) ecorr3_nucl=energia(38) + ecation_prot=energia(41) + ecationcation=energia(42) + escbase=energia(46) + epepbase=energia(47) + escpho=energia(48) + epeppho=energia(49) +! energia(41)=ecation_prot +! energia(42)=ecationcation #ifdef SPLITELE @@ -787,7 +849,9 @@ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl& +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb& +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl& - +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_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 #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc & @@ -799,7 +863,9 @@ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl& +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb& +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl& - +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_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 #endif energia(0)=etot ! detecting NaNQ @@ -927,6 +993,8 @@ real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,& ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,& ecorr3_nucl + real(kind=8) :: ecation_prot,ecationcation + real(kind=8) :: escbase,epepbase,escpho,epeppho etot=energia(0) evdw=energia(1) @@ -973,7 +1041,12 @@ etors_d_nucl=energia(36) ecorr_nucl=energia(37) ecorr3_nucl=energia(38) - + ecation_prot=energia(41) + ecationcation=energia(42) + escbase=energia(46) + epepbase=energia(47) + escpho=energia(48) + epeppho=energia(49) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,& estr,wbond,ebe,wang,& @@ -987,7 +1060,8 @@ evdwpp,wvdwpp_nucl,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb,& evdwsb,wvdwsb,eelsb,welsb,esbloc,wsbloc,etors_nucl,wtor_nucl,& etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,& - ecorr3_nucl,wcorr3_nucl, & + ecorr3_nucl,wcorr3_nucl,ecation_prot,wcatprot,ecationcation,wcatcat, & + escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,& etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & @@ -1029,6 +1103,12 @@ 'ETORSD_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(double torsional)'/ & 'ECORR_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(multibody 4th order)'/ & 'ECORR3_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(multibody 3th order)'/ & + 'ECATPROT=',1pE16.6,' WEIGHT=',1pD16.6,'(cation prot)'/ & + 'ECATCAT=',1pE16.6,' WEIGHT=',1pD16.6,'(cation cation)'/ & + 'ESCBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(sc-prot nucl-base)'/ & + 'EPEPBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-base)'/ & + 'ESCPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(sc-prot nucl-phosphate)'/& + 'EPEPPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-phosphate)'/& 'ETOT= ',1pE16.6,' (total)') #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,& @@ -1043,7 +1123,8 @@ evdwpp,wvdwpp_nucl,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb& evdwsb,wvdwsb,eelsb,welsb,esbloc,wsbloc,etors_nucl,wtor_nucl& etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,& - ecorr3_nucl,wcorr3_nucl, & + ecorr3_nucl,wcorr3_nucl,ecation_prot,wcatprot,ecationcation,wcatcat, & + escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,& etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & @@ -1084,6 +1165,12 @@ 'ETORSD_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(double torsional)'/ & 'ECORR_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(multibody 4th order)'/ & 'ECORR3_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(multibody 3th order)'/ & + 'ECATPROT=',1pE16.6,' WEIGHT=',1pD16.6,'(cation prot)'/ & + 'ECATCAT=',1pE16.6,' WEIGHT=',1pD16.6,'(cation cation)'/ & + 'ESCBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(sc-prot nucl-base)'/ & + 'EPEPBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-base)'/ & + 'ESCPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(sc-prot nucl-phosphate)'/& + 'EPEPPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-phosphate)'/& 'ETOT= ',1pE16.6,' (total)') #endif return @@ -2492,13 +2579,21 @@ endif ! if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then if (i.gt. nnt+2 .and. i.lt.nct+2) then + if (itype(i-2,1).eq.0) then + iti=ntortyp+1 + else iti = itortyp(itype(i-2,1)) + endif else iti=ntortyp+1 endif ! if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then + if (itype(i-1,1).eq.0) then + iti1=ntortyp+1 + else iti1 = itortyp(itype(i-1,1)) + endif else iti1=ntortyp+1 endif @@ -2539,7 +2634,9 @@ enddo ! if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then - if (itype(i-1,1).le.ntyp) then + if (itype(i-1,1).eq.0) then + iti1=ntortyp+1 + elseif (itype(i-1,1).le.ntyp) then iti1 = itortyp(itype(i-1,1)) else iti1=ntortyp+1 @@ -6690,12 +6787,12 @@ difi=phii-phi0(i) if (difi.gt.drange(i)) then difi=difi-drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 else if (difi.lt.-drange(i)) then difi=difi+drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 endif ! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii, ! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) @@ -6801,12 +6898,12 @@ difi=pinorm(phii-phi0(i)) if (difi.gt.drange(i)) then difi=difi-drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 else if (difi.lt.-drange(i)) then difi=difi+drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 else difi=0.0 endif @@ -10448,7 +10545,17 @@ wvdwpsb*(gvdwpsb(j,i)+gvdwpsb1(j,i))+ & wvdwsb*gvdwsbc(j,i)+welsb*gelsbc(j,i)+ & wcorr_nucl*gradcorr_nucl(j,i)& - +wcorr3_nucl*gradcorr3_nucl(j,i) + +wcorr3_nucl*gradcorr3_nucl(j,i)+& + wcatprot* gradpepcat(j,i)+ & + wcatcat*gradcatcat(j,i)+ & + wscbase*gvdwc_scbase(j,i)+ & + wpepbase*gvdwc_pepbase(j,i)+& + wscpho*gvdwc_scpho(j,i)+ & + wpeppho*gvdwc_peppho(j,i) + + + + enddo enddo @@ -10475,8 +10582,16 @@ +wvdwpp_nucl*gvdwpp_nucl(j,i)+welpp*gelpp(j,i)+ & wvdwpsb*(gvdwpsb(j,i)+gvdwpsb1(j,i))+ & wvdwsb*gvdwsbc(j,i)+welsb*gelsbc(j,i)+ & - wcorr_nucl*gradcorr_nucl(j,i) - +wcorr3_nucl*gradcorr3_nucl(j,i) + wcorr_nucl*gradcorr_nucl(j,i) & + +wcorr3_nucl*gradcorr3_nucl(j,i) +& + wcatprot* gradpepcat(j,i)+ & + wcatcat*gradcatcat(j,i)+ & + wscbase*gvdwc_scbase(j,i) & + wpepbase*gvdwc_pepbase(j,i)+& + wscpho*gvdwc_scpho(j,i)+& + wpeppho*gvdwc_peppho(j,i) + + enddo enddo #endif @@ -10734,7 +10849,13 @@ +welsb*gelsbx(j,i) & +wcorr_nucl*gradxorr_nucl(j,i)& +wcorr3_nucl*gradxorr3_nucl(j,i) & - +wsbloc*gsblocx(j,i) + +wsbloc*gsblocx(j,i) & + +wcatprot* gradpepcatx(j,i)& + +wscbase*gvdwx_scbase(j,i) & + +wpepbase*gvdwx_pepbase(j,i)& + +wscpho*gvdwx_scpho(j,i) +! if (i.eq.3) print *,"tu?", wscpho,gvdwx_scpho(j,i) + enddo enddo #ifdef DEBUG @@ -16169,9 +16290,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' #ifdef TIMING time01=MPI_Wtime() #endif - print *,"gcart_two",gxcart(2,2),gradx(2,2,icg) +! print *,"gcart_two",gxcart(2,2),gradx(2,2,icg) call int_to_cart - print *,"gcart_two",gxcart(2,2),gradx(2,2,icg) +! print *,"gcart_two",gxcart(2,2),gradx(2,2,icg) #ifdef TIMING time_inttocart=time_inttocart+MPI_Wtime()-time01 @@ -16338,6 +16459,16 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' gelsbx(j,i)=0.0d0 gsbloc(j,i)=0.0d0 gsblocx(j,i)=0.0d0 + gradpepcat(j,i)=0.0d0 + gradpepcatx(j,i)=0.0d0 + gradcatcat(j,i)=0.0d0 + gvdwx_scbase(j,i)=0.0d0 + gvdwc_scbase(j,i)=0.0d0 + gvdwx_pepbase(j,i)=0.0d0 + gvdwc_pepbase(j,i)=0.0d0 + gvdwx_scpho(j,i)=0.0d0 + gvdwc_scpho(j,i)=0.0d0 + gvdwc_peppho(j,i)=0.0d0 enddo enddo do i=0,nres @@ -19717,7 +19848,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' allocate(gradcorr3_nucl(3,-1:nres)) allocate(gradxorr3_nucl(3,-1:nres)) allocate(gvdwpp_nucl(3,-1:nres)) - + allocate(gradpepcat(3,-1:nres)) + allocate(gradpepcatx(3,-1:nres)) + allocate(gradcatcat(3,-1:nres)) !(3,maxres) allocate(grad_shield_side(3,50,nres)) allocate(grad_shield_loc(3,50,nres)) @@ -19746,6 +19879,14 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !(3,maxres) allocate(gsccor_loc(-1:nres)) !(maxres) + allocate(gvdwx_scbase(3,-1:nres)) + allocate(gvdwc_scbase(3,-1:nres)) + allocate(gvdwx_pepbase(3,-1:nres)) + allocate(gvdwc_pepbase(3,-1:nres)) + allocate(gvdwx_scpho(3,-1:nres)) + allocate(gvdwc_scpho(3,-1:nres)) + allocate(gvdwc_peppho(3,-1:nres)) + allocate(dtheta(3,2,-1:nres)) !(3,2,maxres) allocate(gscloc(3,-1:nres)) @@ -20666,8 +20807,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' evdwij=evdwij*eps2rt evdwsb=evdwsb+evdwij if (lprn) then - sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) - epsi=bb(itypi,itypj)**2/aa(itypi,itypj) + sigm=dabs(aa_nucl(itypi,itypj)/bb_nucl(itypi,itypj))**(1.0D0/6.0D0) + epsi=bb_nucl(itypi,itypj)**2/aa_nucl(itypi,itypj) write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi,2),i,restyp(itypj,2),j, & epsi,sigm,chi1,chi2,chip1,chip2, & @@ -21622,15 +21763,19 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' epscalc=0.05 r06 = rcat0**6 r012 = r06**2 - k0 = 332*(2*2)/80 + k0 = 332.0*(2.0*2.0)/80.0 itmp=0 + do i=1,4 itmp=itmp+nres_molec(i) enddo - do i=itmp+1,itmp+nres_molec(i)-1 +! write(iout,*) "itmp",itmp + do i=itmp+1,itmp+nres_molec(5)-1 + xi=c(1,i) yi=c(2,i) zi=c(3,i) + xi=mod(xi,boxxsize) if (xi.lt.0) xi=xi+boxxsize yi=mod(yi,boxysize) @@ -21638,7 +21783,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' zi=mod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize - do j=i+1,itmp+nres_molec(i) + do j=i+1,itmp+nres_molec(5) +! print *,i,j,'catcat' xj=c(1,j) yj=c(2,j) zj=c(3,j) @@ -21648,6 +21794,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' if (yj.lt.0) yj=yj+boxysize zj=dmod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize +! write(iout,*) c(1,i),xi,xj,"xy",boxxsize dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 xj_safe=xj yj_safe=yj @@ -21701,10 +21848,11 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' enddo do k=1,3 gg(k) = dEvan1Cmcat(k)+dEvan2Cmcat(k)+dEeleccat(k) - gradcatcat(k,i)=gradcatcat(k,i)+gg(k) - gradcatcat(k,j)=gradcatcat(k,j)-gg(k) + gradcatcat(k,i)=gradcatcat(k,i)-gg(k) + gradcatcat(k,j)=gradcatcat(k,j)+gg(k) enddo +! write(iout,*) "ecatcat",i,j, ecationcation,xj,yj,zj ecationcation=ecationcation+Evan1cat+Evan2cat+Eeleccat enddo enddo @@ -21712,20 +21860,51 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' end subroutine ecatcat !--------------------------------------------------------------------------- subroutine ecat_prot(ecation_prot) - integer i,j,k,subchap,itmp + 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 + 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 real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,& - gg,r + 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 + wconst=78 + 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 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=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)) @@ -21778,14 +21957,2149 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 + 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) + if (yi.lt.0) yi=yi+boxysize + zi=mod(zi,boxzsize) + if (zi.lt.0) zi=zi+boxzsize + do k=1,3 + cm1(k)=dc(k,i+nres) + enddo + cm1mag=sqrt(cm1(1)**2+cm1(2)**2+cm1(3)**2) + do j=itmp+1,itmp+nres_molec(5) + 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 + if(itype(i,1).eq.15.or.itype(i,1).eq.16) then + if(itype(i,1).eq.16) then + inum=1 + else + inum=2 + endif + do k=1,6 + vcatprm(k)=catprm(k,inum) + enddo + dASGL=catprm(7,inum) + do k=1,3 + vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres) + valpha(k)=c(k,i) + vcat(k)=c(k,j) + enddo + do k=1,3 + dx(k) = vcat(k)-vcm(k) + enddo + do k=1,3 + v1(k)=(vcm(k)-valpha(k)) + v2(k)=(vcat(k)-valpha(k)) + enddo + v1m = sqrt(v1(1)**2+v1(2)**2+v1(3)**2) + v2m = sqrt(v2(1)**2+v2(2)**2+v2(3)**2) + v1dpv2 = v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3) + +! The weights of the energy function calculated from +!The quantum mechanical GAMESS simulations of calcium with ASP/GLU + wh2o=78 + wc = vcatprm(1) + wc=wc/wh2o + wdip =vcatprm(2) + wdip=wdip/wh2o + wquad1 =vcatprm(3) + wquad1=wquad1/wh2o + wquad2 = vcatprm(4) + wquad2=wquad2/wh2o + wquad2p = 1-wquad2 + wvan1 = vcatprm(5) + wvan2 =vcatprm(6) + opt = dx(1)**2+dx(2)**2 + rsecp = opt+dx(3)**2 + rs = sqrt(rsecp) + rthrp = rsecp*rs + rfourp = rthrp*rs + rsixp = rfourp*rsecp + reight=rsixp*rsecp + Ir = 1.0d0/rs + Irsecp = 1/rsecp + Irthrp = Irsecp/rs + Irfourp = Irthrp/rs + Irsixp = 1/rsixp + Ireight=1/reight + Irtw=Irsixp*Irsixp + Irthir=Irtw/rs + Irfourt=Irthir/rs + opt1 = (4*rs*dx(3)*wdip) + opt2 = 6*rsecp*wquad1*opt + opt3 = wquad1*wquad2p*Irsixp + opt4 = (wvan1*wvan2**12) + opt5 = opt4*12*Irfourt + opt6 = 2*wvan1*wvan2**6 + opt7 = 6*opt6*Ireight + opt8 = wdip/v1m + opt10 = wdip/v2m + opt11 = (rsecp*v2m)**2 + opt12 = (rsecp*v1m)**2 + opt14 = (v1m*v2m*rsecp)**2 + opt15 = -wquad1/v2m**2 + opt16 = (rthrp*(v1m*v2m)**2)**2 + opt17 = (v1m**2*rthrp)**2 + opt18 = -wquad1/rthrp + opt19 = (v1m**2*v2m**2)**2 + Ec = wc*Ir + do k=1,3 + dEcCat(k) = -(dx(k)*wc)*Irthrp + dEcCm(k)=(dx(k)*wc)*Irthrp + dEcCalp(k)=0.0d0 + enddo + Edip=opt8*(v1dpv2)/(rsecp*v2m) + do k=1,3 + dEdipCat(k)=opt8*(v1(k)*rsecp*v2m-((v2(k)/v2m & + *rsecp+2*dx(k)*v2m)*v1dpv2))/opt11 + dEdipCm(k)=opt10*(v2(k)*rsecp*v1m-((v1(k)/v1m & + *rsecp-2*dx(k)*v1m)*v1dpv2))/opt12 + dEdipCalp(k)=wdip*((-v1(k)-v2(k))*rsecp*v1m & + *v2m-(-v1(k)/v1m*v2m*rsecp-v2(k)/v2m*v1m*rsecp) & + *v1dpv2)/opt14 + enddo + Equad1=-wquad1*v1dpv2**2/(rthrp*(v1m*v2m)**2) + do k=1,3 + dEquad1Cat(k)=-wquad1*(2*v1(k)*v1dpv2*(rthrp* & + (v1m*v2m)**2)-(3*dx(k)*rs*(v1m*v2m)**2+2*v1m*2* & + v2(k)*1/2*1/v2m*v1m*v2m*rthrp)*v1dpv2**2)/opt16 + dEquad1Cm(k)=-wquad1*(2*v2(k)*v1dpv2*(rthrp* & + (v1m*v2m)**2)-(-3*dx(k)*rs*(v1m*v2m)**2+2*v2m*2* & + v1(k)*1/2*1/v1m*v2m*v1m*rthrp)*v1dpv2**2)/opt16 + dEquad1Calp(k)=opt18*(2*(-v1(k)-v2(k))*v1dpv2* & + v1m**2*v2m**2-(-2*v1(k)*v2m**2-2*v2(k)*v1m**2)* & + v1dpv2**2)/opt19 + enddo + Equad2=wquad1*wquad2p*Irthrp + do k=1,3 + dEquad2Cat(k)=-3*dx(k)*rs*opt3 + dEquad2Cm(k)=3*dx(k)*rs*opt3 + dEquad2Calp(k)=0.0d0 + enddo + Evan1=opt4*Irtw + do k=1,3 + dEvan1Cat(k)=-dx(k)*opt5 + dEvan1Cm(k)=dx(k)*opt5 + dEvan1Calp(k)=0.0d0 + enddo + Evan2=-opt6*Irsixp + do k=1,3 + dEvan2Cat(k)=dx(k)*opt7 + dEvan2Cm(k)=-dx(k)*opt7 + dEvan2Calp(k)=0.0d0 + enddo + ecation_prot=ecation_prot+Ec+Edip+Equad1+Equad2+Evan1+Evan2 +! print *,ecation_prot,Ec+Edip+Equad1+Equad2+Evan1+Evan2 + + do k=1,3 + dEtotalCat(k)=dEcCat(k)+dEdipCat(k)+dEquad1Cat(k)+ & + dEquad2Cat(k)+dEvan1Cat(k)+dEvan2Cat(k) +!c write(*,*) 'dEtotalCat inside', (dEtotalCat(l),l=1,3) + dEtotalCm(k)=dEcCm(k)+dEdipCm(k)+dEquad1Cm(k)+ & + dEquad2Cm(k)+dEvan1Cm(k)+dEvan2Cm(k) + dEtotalCalp(k)=dEcCalp(k)+dEdipCalp(k)+dEquad1Calp(k) & + +dEquad2Calp(k)+dEvan1Calp(k)+dEvan2Calp(k) + enddo + dscmag = 0.0d0 + do k=1,3 + dscvec(k) = dc(k,i+nres) + dscmag = dscmag+dscvec(k)*dscvec(k) + enddo + dscmag3 = dscmag + dscmag = sqrt(dscmag) + dscmag3 = dscmag3*dscmag + constA = 1.0d0+dASGL/dscmag + constB = 0.0d0 + do k=1,3 + constB = constB+dscvec(k)*dEtotalCm(k) + enddo + constB = constB*dASGL/dscmag3 + do k=1,3 + gg(k) = dEtotalCm(k)+dEtotalCalp(k) + gradpepcatx(k,i)=gradpepcatx(k,i)+ & + constA*dEtotalCm(k)-constB*dscvec(k) +! print *,j,constA,dEtotalCm(k),constB,dscvec(k) + gradpepcat(k,i)=gradpepcat(k,i)+gg(k) + gradpepcat(k,j)=gradpepcat(k,j)+dEtotalCat(k) + enddo + else if (itype(i,1).eq.13.or.itype(i,1).eq.14) then + if(itype(i,1).eq.14) then + inum=3 + else + inum=4 + endif + do k=1,6 + vcatprm(k)=catprm(k,inum) + enddo + dASGL=catprm(7,inum) + do k=1,3 + vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres) + valpha(k)=c(k,i) + vcat(k)=c(k,j) + enddo - + do k=1,3 + dx(k) = vcat(k)-vcm(k) + enddo + do k=1,3 + v1(k)=(vcm(k)-valpha(k)) + v2(k)=(vcat(k)-valpha(k)) + enddo + v1m = sqrt(v1(1)**2+v1(2)**2+v1(3)**2) + v2m = sqrt(v2(1)**2+v2(2)**2+v2(3)**2) + v1dpv2 = v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3) +! The weights of the energy function calculated from +!The quantum mechanical GAMESS simulations of ASN/GLN with calcium + wh2o=78 + wdip =vcatprm(2) + wdip=wdip/wh2o + wquad1 =vcatprm(3) + wquad1=wquad1/wh2o + wquad2 = vcatprm(4) + wquad2=wquad2/wh2o + wquad2p = 1-wquad2 + wvan1 = vcatprm(5) + wvan2 =vcatprm(6) + opt = dx(1)**2+dx(2)**2 + rsecp = opt+dx(3)**2 + rs = sqrt(rsecp) + rthrp = rsecp*rs + rfourp = rthrp*rs + rsixp = rfourp*rsecp + reight=rsixp*rsecp + Ir = 1.0d0/rs + Irsecp = 1/rsecp + Irthrp = Irsecp/rs + Irfourp = Irthrp/rs + Irsixp = 1/rsixp + Ireight=1/reight + Irtw=Irsixp*Irsixp + Irthir=Irtw/rs + Irfourt=Irthir/rs + opt1 = (4*rs*dx(3)*wdip) + opt2 = 6*rsecp*wquad1*opt + opt3 = wquad1*wquad2p*Irsixp + opt4 = (wvan1*wvan2**12) + opt5 = opt4*12*Irfourt + opt6 = 2*wvan1*wvan2**6 + opt7 = 6*opt6*Ireight + opt8 = wdip/v1m + opt10 = wdip/v2m + opt11 = (rsecp*v2m)**2 + opt12 = (rsecp*v1m)**2 + opt14 = (v1m*v2m*rsecp)**2 + opt15 = -wquad1/v2m**2 + opt16 = (rthrp*(v1m*v2m)**2)**2 + opt17 = (v1m**2*rthrp)**2 + opt18 = -wquad1/rthrp + opt19 = (v1m**2*v2m**2)**2 + Edip=opt8*(v1dpv2)/(rsecp*v2m) + do k=1,3 + dEdipCat(k)=opt8*(v1(k)*rsecp*v2m-((v2(k)/v2m& + *rsecp+2*dx(k)*v2m)*v1dpv2))/opt11 + dEdipCm(k)=opt10*(v2(k)*rsecp*v1m-((v1(k)/v1m& + *rsecp-2*dx(k)*v1m)*v1dpv2))/opt12 + dEdipCalp(k)=wdip*((-v1(k)-v2(k))*rsecp*v1m& + *v2m-(-v1(k)/v1m*v2m*rsecp-v2(k)/v2m*v1m*rsecp)& + *v1dpv2)/opt14 + enddo + Equad1=-wquad1*v1dpv2**2/(rthrp*(v1m*v2m)**2) + do k=1,3 + dEquad1Cat(k)=-wquad1*(2*v1(k)*v1dpv2*(rthrp*& + (v1m*v2m)**2)-(3*dx(k)*rs*(v1m*v2m)**2+2*v1m*2*& + v2(k)*1/2*1/v2m*v1m*v2m*rthrp)*v1dpv2**2)/opt16 + dEquad1Cm(k)=-wquad1*(2*v2(k)*v1dpv2*(rthrp*& + (v1m*v2m)**2)-(-3*dx(k)*rs*(v1m*v2m)**2+2*v2m*2*& + v1(k)*1/2*1/v1m*v2m*v1m*rthrp)*v1dpv2**2)/opt16 + dEquad1Calp(k)=opt18*(2*(-v1(k)-v2(k))*v1dpv2* & + v1m**2*v2m**2-(-2*v1(k)*v2m**2-2*v2(k)*v1m**2)*& + v1dpv2**2)/opt19 + enddo + Equad2=wquad1*wquad2p*Irthrp + do k=1,3 + dEquad2Cat(k)=-3*dx(k)*rs*opt3 + dEquad2Cm(k)=3*dx(k)*rs*opt3 + dEquad2Calp(k)=0.0d0 + enddo + Evan1=opt4*Irtw + do k=1,3 + dEvan1Cat(k)=-dx(k)*opt5 + dEvan1Cm(k)=dx(k)*opt5 + dEvan1Calp(k)=0.0d0 + enddo + Evan2=-opt6*Irsixp + do k=1,3 + dEvan2Cat(k)=dx(k)*opt7 + dEvan2Cm(k)=-dx(k)*opt7 + dEvan2Calp(k)=0.0d0 + enddo + ecation_prot = ecation_prot+Edip+Equad1+Equad2+Evan1+Evan2 + do k=1,3 + dEtotalCat(k)=dEdipCat(k)+dEquad1Cat(k)+ & + dEquad2Cat(k)+dEvan1Cat(k)+dEvan2Cat(k) + dEtotalCm(k)=dEdipCm(k)+dEquad1Cm(k)+ & + dEquad2Cm(k)+dEvan1Cm(k)+dEvan2Cm(k) + dEtotalCalp(k)=dEdipCalp(k)+dEquad1Calp(k) & + +dEquad2Calp(k)+dEvan1Calp(k)+dEvan2Calp(k) + enddo + dscmag = 0.0d0 + do k=1,3 + dscvec(k) = c(k,i+nres)-c(k,i) + dscmag = dscmag+dscvec(k)*dscvec(k) + enddo + dscmag3 = dscmag + dscmag = sqrt(dscmag) + dscmag3 = dscmag3*dscmag + constA = 1+dASGL/dscmag + constB = 0.0d0 + do k=1,3 + constB = constB+dscvec(k)*dEtotalCm(k) + enddo + constB = constB*dASGL/dscmag3 + do k=1,3 + gg(k) = dEtotalCm(k)+dEtotalCalp(k) + gradpepcatx(k,i)=gradpepcatx(k,i)+ & + constA*dEtotalCm(k)-constB*dscvec(k) + gradpepcat(k,i)=gradpepcat(k,i)+gg(k) + gradpepcat(k,j)=gradpepcat(k,j)+dEtotalCat(k) + enddo + else + rcal = 0.0d0 + do k=1,3 + r(k) = c(k,j)-c(k,i+nres) + rcal = rcal+r(k)*r(k) + enddo + ract=sqrt(rcal) + rocal=1.5 + epscalc=0.2 + r0p=0.5*(rocal+sig0(itype(i,1))) + r06 = r0p**6 + r012 = r06*r06 + Evan1=epscalc*(r012/rcal**6) + Evan2=epscalc*2*(r06/rcal**3) + r4 = rcal**4 + r7 = rcal**7 + do k=1,3 + dEvan1Cm(k) = 12*r(k)*epscalc*r012/r7 + dEvan2Cm(k) = 12*r(k)*epscalc*r06/r4 + enddo + do k=1,3 + dEtotalCm(k)=dEvan1Cm(k)+dEvan2Cm(k) + enddo + ecation_prot = ecation_prot+ Evan1+Evan2 + do k=1,3 + gradpepcatx(k,i)=gradpepcatx(k,i)+ & + dEtotalCm(k) + gradpepcat(k,i)=gradpepcat(k,i)+dEtotalCm(k) + gradpepcat(k,j)=gradpepcat(k,j)-dEtotalCm(k) + enddo + endif ! 13-16 residues + enddo !j + enddo !i return end subroutine ecat_prot !---------------------------------------------------------------------------- !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- + subroutine eprot_sc_base(escbase) + use calc_data +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +! include 'COMMON.GEO' +! include 'COMMON.VAR' +! include 'COMMON.LOCAL' +! include 'COMMON.CHAIN' +! include 'COMMON.DERIV' +! include 'COMMON.NAMES' +! include 'COMMON.INTERACT' +! include 'COMMON.IOUNITS' +! include 'COMMON.CALC' +! include 'COMMON.CONTROL' +! include 'COMMON.SBRIDGE' + logical :: lprn +!el local variables + integer :: iint,itypi,itypi1,itypj,subchap + real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi + real(kind=8) :: evdw,sig0ij + real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,& + dist_temp, dist_init,aa,bb,ssgradlipi,ssgradlipj, & + sslipi,sslipj,faclip + integer :: ii + real(kind=8) :: fracinbuf + real (kind=8) :: escbase + real (kind=8),dimension(4):: ener + real(kind=8) :: b1,b2,b3,b4,egb,eps_in,eps_inout_fac,eps_out + real(kind=8) :: ECL,Elj,Equad,Epol,eheadtail,rhead,dGCLOM2,& + sqom1,sqom2,sqom12,c1,c2,c3,pom,Lambf,sparrow,& + Chif,ChiLambf,bat,eagle,top,bot,botsq,Fcav,dtop,dFdR,dFdOM1,& + dFdOM2,w1,w2,w3,dGCLdR,dFdL,dFdOM12,dbot ,& + r1,eps_head,alphapol1,pis,facd2,d2,facd1,d1,erdxj,erdxi,federmaus,& + dPOLdR1,dFGBdOM2,dFGBdR1,dPOLdFGB1,RR1,MomoFac1,hawk,d1i,d1j,& + sig1,sig2,chis12,chis2,ee1,fgb1,a12sq,chis1 + real(kind=8),dimension(3,2)::chead,erhead_tail + real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead + integer troll + eps_out=80.0d0 + escbase=0.0d0 +! do i=1,nres_molec(1) + do i=ibond_start,ibond_end + if (itype(i,1).eq.ntyp1_molec(1)) cycle + itypi = itype(i,1) + dxi = dc_norm(1,nres+i) + dyi = dc_norm(2,nres+i) + dzi = dc_norm(3,nres+i) + dsci_inv = vbld_inv(i+nres) + xi=c(1,nres+i) + yi=c(2,nres+i) + zi=c(3,nres+i) + 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=nres_molec(1)+1,nres_molec(2)+nres_molec(1) + itypj= itype(j,2) + if (itype(j,2).eq.ntyp1_molec(2))cycle + xj=c(1,j+nres) + yj=c(2,j+nres) + zj=c(3,j+nres) + 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 + dxj = dc_norm( 1, nres+j ) + dyj = dc_norm( 2, nres+j ) + dzj = dc_norm( 3, nres+j ) +! print *,i,j,itypi,itypj + d1i = dhead_scbasei(itypi,itypj) !this is shift of dipole/charge + d1j = dhead_scbasej(itypi,itypj) !this is shift of dipole/charge +! d1i=0.0d0 +! d1j=0.0d0 +! BetaT = 1.0d0 / (298.0d0 * Rb) +! Gay-berne var's + sig0ij = sigma_scbase( itypi,itypj ) + chi1 = chi_scbase( itypi, itypj,1 ) + chi2 = chi_scbase( itypi, itypj,2 ) +! chi1=0.0d0 +! chi2=0.0d0 + chi12 = chi1 * chi2 + chip1 = chipp_scbase( itypi, itypj,1 ) + chip2 = chipp_scbase( itypi, itypj,2 ) +! chip1=0.0d0 +! chip2=0.0d0 + chip12 = chip1 * chip2 +! not used by momo potential, but needed by sc_angular which is shared +! by all energy_potential subroutines + alf1 = 0.0d0 + alf2 = 0.0d0 + alf12 = 0.0d0 + a12sq = rborn_scbasei(itypi,itypj) * rborn_scbasej(itypi,itypj) +! a12sq = a12sq * a12sq +! charge of amino acid itypi is... + chis1 = chis_scbase(itypi,itypj,1) + chis2 = chis_scbase(itypi,itypj,2) + chis12 = chis1 * chis2 + sig1 = sigmap1_scbase(itypi,itypj) + sig2 = sigmap2_scbase(itypi,itypj) +! write (*,*) "sig1 = ", sig1 +! write (*,*) "sig2 = ", sig2 +! alpha factors from Fcav/Gcav + b1 = alphasur_scbase(1,itypi,itypj) +! b1=0.0d0 + b2 = alphasur_scbase(2,itypi,itypj) + b3 = alphasur_scbase(3,itypi,itypj) + b4 = alphasur_scbase(4,itypi,itypj) +! used to determine whether we want to do quadrupole calculations +! used by Fgb + eps_in = epsintab_scbase(itypi,itypj) + if (eps_in.eq.0.0) eps_in=1.0 + eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out)) +! write (*,*) "eps_inout_fac = ", eps_inout_fac +!------------------------------------------------------------------- +! tail location and distance calculations + 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+nres) + d1i * dc_norm(k, i+nres) + chead(k,2) = c(k, j+nres) + d1j * dc_norm(k, j+nres) +! 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 = 1.0/rij - 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_scbase(itypi,itypj) +! c1 = 0.0d0 + c2 = fac * bb_scbase(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 + c1 = c1 * eps1 * eps2rt**2 * eps3rt**2 + fac = -expon * (c1 + evdwij) * rij_shift + sigder = fac * sigder +! fac = rij * fac +! Calculate distance derivative + gg(1) = fac + gg(2) = fac + gg(3) = fac +! if (b2.gt.0.0) then + fac = chis1 * sqom1 + chis2 * sqom2 & + - 2.0d0 * chis12 * om1 * om2 * om12 +! we will use pom later in Gcav, so dont mess with it! + pom = 1.0d0 - chis1 * chis2 * sqom12 + Lambf = (1.0d0 - (fac / pom)) + Lambf = dsqrt(Lambf) + sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0) +! write (*,*) "sparrow = ", sparrow + Chif = 1.0d0/rij * sparrow + ChiLambf = Chif * Lambf + eagle = dsqrt(ChiLambf) + bat = ChiLambf ** 11.0d0 + top = b1 * ( eagle + b2 * ChiLambf - b3 ) + bot = 1.0d0 + b4 * (ChiLambf ** 12.0d0) + botsq = bot * bot + Fcav = top / bot +! print *,i,j,Fcav + dtop = b1 * ((Lambf / (2.0d0 * eagle)) + (b2 * Lambf)) + dbot = 12.0d0 * b4 * bat * Lambf + dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow +! dFdR = 0.0d0 +! write (*,*) "dFcav/dR = ", dFdR + dtop = b1 * ((Chif / (2.0d0 * eagle)) + (b2 * Chif)) + dbot = 12.0d0 * b4 * 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) +! dFdL = 0.0d0 + dCAVdOM1 = dFdL * ( dFdOM1 ) + dCAVdOM2 = dFdL * ( dFdOM2 ) + dCAVdOM12 = dFdL * ( dFdOM12 ) + + ertail(1) = xj*rij + ertail(2) = yj*rij + ertail(3) = zj*rij +! eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 +! eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 +! eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 & +! -2.0D0*alf12*eps3der+sigder*sigsq_om12 +! print *,"EOMY",eom1,eom2,eom12 +! erdxi = scalar( ertail(1), dC_norm(1,i+nres) ) +! erdxj = scalar( ertail(1), dC_norm(1,j+nres) ) +! here dtail=0.0 +! facd1 = dtail(1,itypi,itypj) * vbld_inv(i+nres) +! facd2 = dtail(2,itypi,itypj) * vbld_inv(j+nres) + DO k = 1, 3 +! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i) +! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j) + pom = ertail(k) +!-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres)) + gvdwx_scbase(k,i) = gvdwx_scbase(k,i) & + - (( dFdR + gg(k) ) * pom) +! +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) & +! +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv +! & - ( dFdR * pom ) + pom = ertail(k) +!-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres)) + gvdwx_scbase(k,j) = gvdwx_scbase(k,j) & + + (( dFdR + gg(k) ) * pom) +! +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) & +! +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv +!c! & + ( dFdR * pom ) + + gvdwc_scbase(k,i) = gvdwc_scbase(k,i) & + - (( dFdR + gg(k) ) * ertail(k)) +!c! & - ( dFdR * ertail(k)) + + gvdwc_scbase(k,j) = gvdwc_scbase(k,j) & + + (( dFdR + gg(k) ) * ertail(k)) +!c! & + ( dFdR * ertail(k)) + + gg(k) = 0.0d0 +!c! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i) +!c! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j) + END DO + +! else + +! endif +!Now dipole-dipole + if (wdipdip_scbase(2,itypi,itypj).gt.0.0d0) then + w1 = wdipdip_scbase(1,itypi,itypj) + w2 = -wdipdip_scbase(3,itypi,itypj)/2.0 + w3 = wdipdip_scbase(2,itypi,itypj) +!c!------------------------------------------------------------------- +!c! ECL + fac = (om12 - 3.0d0 * om1 * om2) + c1 = (w1 / (Rhead**3.0d0)) * fac + c2 = (w2 / Rhead ** 6.0d0) & + * (4.0d0 + fac * fac -3.0d0 * (sqom1 + sqom2)) + c3= (w3/ Rhead ** 6.0d0) & + * (2.0d0 - 2.0d0*fac*fac +3.0d0*(sqom1 + sqom2)) + ECL = c1 - c2 + c3 +!c! write (*,*) "w1 = ", w1 +!c! write (*,*) "w2 = ", w2 +!c! write (*,*) "om1 = ", om1 +!c! write (*,*) "om2 = ", om2 +!c! write (*,*) "om12 = ", om12 +!c! write (*,*) "fac = ", fac +!c! write (*,*) "c1 = ", c1 +!c! write (*,*) "c2 = ", c2 +!c! write (*,*) "Ecl = ", Ecl +!c! write (*,*) "c2_1 = ", (w2 / Rhead ** 6.0d0) +!c! write (*,*) "c2_2 = ", +!c! & (4.0d0 + fac * fac -3.0d0 * (sqom1 + sqom2)) +!c!------------------------------------------------------------------- +!c! dervative of ECL is GCL... +!c! dECL/dr + c1 = (-3.0d0 * w1 * fac) / (Rhead ** 4.0d0) + c2 = (-6.0d0 * w2) / (Rhead ** 7.0d0) & + * (4.0d0 + fac * fac - 3.0d0 * (sqom1 + sqom2)) + c3= (-6.0d0 * w3) / (Rhead ** 7.0d0) & + * (2.0d0 - 2.0d0*fac*fac +3.0d0*(sqom1 + sqom2)) + dGCLdR = c1 - c2 + c3 +!c! dECL/dom1 + c1 = (-3.0d0 * w1 * om2 ) / (Rhead**3.0d0) + c2 = (-6.0d0 * w2) / (Rhead**6.0d0) & + * ( om2 * om12 - 3.0d0 * om1 * sqom2 + om1 ) + c3 =(6.0d0*w3/ Rhead ** 6.0d0)*(om1-2.0d0*(fac)*(-om2)) + dGCLdOM1 = c1 - c2 + c3 +!c! dECL/dom2 + c1 = (-3.0d0 * w1 * om1 ) / (Rhead**3.0d0) + c2 = (-6.0d0 * w2) / (Rhead**6.0d0) & + * ( om1 * om12 - 3.0d0 * sqom1 * om2 + om2 ) + c3 =(6.0d0*w3/ Rhead ** 6.0d0)*(om2-2.0d0*(fac)*(-om1)) + dGCLdOM2 = c1 - c2 + c3 +!c! dECL/dom12 + c1 = w1 / (Rhead ** 3.0d0) + c2 = ( 2.0d0 * w2 * fac ) / Rhead ** 6.0d0 + c3 = (w3/ Rhead ** 6.0d0)*(-4.0d0*fac) + dGCLdOM12 = c1 - c2 + c3 + DO k= 1, 3 + erhead(k) = Rhead_distance(k)/Rhead + END DO + erdxi = scalar( erhead(1), dC_norm(1,i+nres) ) + erdxj = scalar( erhead(1), dC_norm(1,j+nres) ) + facd1 = d1i * vbld_inv(i+nres) + facd2 = d1j * vbld_inv(j+nres) + DO k = 1, 3 + + pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres)) + gvdwx_scbase(k,i) = gvdwx_scbase(k,i) & + - dGCLdR * pom + pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres)) + gvdwx_scbase(k,j) = gvdwx_scbase(k,j) & + + dGCLdR * pom + + gvdwc_scbase(k,i) = gvdwc_scbase(k,i) & + - dGCLdR * erhead(k) + gvdwc_scbase(k,j) = gvdwc_scbase(k,j) & + + dGCLdR * erhead(k) + END DO + endif +!now charge with dipole eg. ARG-dG + if (wqdip_scbase(2,itypi,itypj).gt.0.0d0) then + alphapol1 = alphapol_scbase(itypi,itypj) + w1 = wqdip_scbase(1,itypi,itypj) + w2 = wqdip_scbase(2,itypi,itypj) +! w1=0.0d0 +! w2=0.0d0 +! pis = sig0head_scbase(itypi,itypj) +! eps_head = epshead_scbase(itypi,itypj) +!c!------------------------------------------------------------------- +!c! R1 - distance between head of ith side chain and tail of jth sidechain + R1 = 0.0d0 + DO k = 1, 3 +!c! Calculate head-to-tail distances tail is center of side-chain + R1=R1+(c(k,j+nres)-chead(k,1))**2 + END DO +!c! Pitagoras + R1 = dsqrt(R1) + +!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 * om1 + hawk = w2 * (1.0d0 - sqom2) + 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) / (Rhead**2.0d0) +!c! dF/dom2 + dGCLdOM2 = (2.0d0 * w2 * om2) / (Rhead ** 4.0d0) +!c-------------------------------------------------------------------- +!c Polarization energy +!c Epol + MomoFac1 = (1.0d0 - chi1 * sqom2) + RR1 = R1 * R1 / MomoFac1 + ee1 = exp(-( RR1 / (4.0d0 * a12sq) )) + fgb1 = sqrt( RR1 + a12sq * ee1) +! eps_inout_fac=0.0d0 + epol = 332.0d0 * eps_inout_fac * (( alphapol1 / fgb1 )**4.0d0) +! derivative of Epol is Gpol... + dPOLdFGB1 = -(1328.0d0 * eps_inout_fac * alphapol1 ** 4.0d0) & + / (fgb1 ** 5.0d0) + dFGBdR1 = ( (R1 / MomoFac1) & + * ( 2.0d0 - (0.5d0 * ee1) ) ) & + / ( 2.0d0 * fgb1 ) + dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1)) & + * (2.0d0 - 0.5d0 * ee1) ) & + / (2.0d0 * fgb1) + dPOLdR1 = dPOLdFGB1 * dFGBdR1 +! dPOLdR1 = 0.0d0 + dPOLdOM1 = 0.0d0 + dPOLdOM2 = dPOLdFGB1 * dFGBdOM2 + DO k = 1, 3 + erhead(k) = Rhead_distance(k)/Rhead + erhead_tail(k,1) = ((c(k,j+nres)-chead(k,1))/R1) + END DO + + erdxi = scalar( erhead(1), dC_norm(1,i+nres) ) + erdxj = scalar( erhead(1), dC_norm(1,j+nres) ) + bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) ) +! bat=0.0d0 + federmaus = scalar(erhead_tail(1,1),dC_norm(1,j+nres)) + facd1 = d1i * vbld_inv(i+nres) + facd2 = d1j * vbld_inv(j+nres) +! facd4 = dtail(2,itypi,itypj) * vbld_inv(j+nres) + + DO k = 1, 3 + hawk = (erhead_tail(k,1) + & + facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres))) +! facd1=0.0d0 +! facd2=0.0d0 + pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres)) + gvdwx_scbase(k,i) = gvdwx_scbase(k,i) & + - dGCLdR * pom & + - dPOLdR1 * (erhead_tail(k,1)) +! & - dGLJdR * pom + + pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres)) + gvdwx_scbase(k,j) = gvdwx_scbase(k,j) & + + dGCLdR * pom & + + dPOLdR1 * (erhead_tail(k,1)) +! & + dGLJdR * pom + + + gvdwc_scbase(k,i) = gvdwc_scbase(k,i) & + - dGCLdR * erhead(k) & + - dPOLdR1 * erhead_tail(k,1) +! & - dGLJdR * erhead(k) + + gvdwc_scbase(k,j) = gvdwc_scbase(k,j) & + + dGCLdR * erhead(k) & + + dPOLdR1 * erhead_tail(k,1) +! & + dGLJdR * erhead(k) + + END DO + endif +! print *,i,j,evdwij,epol,Fcav,ECL + escbase=escbase+evdwij+epol+Fcav+ECL + call sc_grad_scbase + enddo + enddo + + return + end subroutine eprot_sc_base + SUBROUTINE sc_grad_scbase + use calc_data + + real (kind=8) :: dcosom1(3),dcosom2(3) + 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 + +! print *,eom1,eom2,eom12,i,j,"eom1,2,12",erij(1),erij(2),erij(3) +! print *,dsci_inv,dscj_inv,dc_norm(2,nres+j),dc_norm(2,nres+i),& +! gg(1),gg(2),"rozne" + DO k = 1, 3 + dcosom1(k) = rij * (dc_norm(k,nres+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) + gvdwx_scbase(k,i)= gvdwx_scbase(k,i) - gg(k) & + + (eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) & + + eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv + gvdwx_scbase(k,j)= gvdwx_scbase(k,j) + gg(k) & + + (eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) & + + eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv + gvdwc_scbase(k,i)=gvdwc_scbase(k,i)-gg(k) + gvdwc_scbase(k,j)=gvdwc_scbase(k,j)+gg(k) + END DO + RETURN + END SUBROUTINE sc_grad_scbase + + + subroutine epep_sc_base(epepbase) + use calc_data + logical :: lprn +!el local variables + integer :: iint,itypi,itypi1,itypj,subchap + real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi + real(kind=8) :: evdw,sig0ij + real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,& + dist_temp, dist_init,aa,bb,ssgradlipi,ssgradlipj, & + sslipi,sslipj,faclip + integer :: ii + real(kind=8) :: fracinbuf + real (kind=8) :: epepbase + real (kind=8),dimension(4):: ener + real(kind=8) :: b1,b2,b3,b4,egb,eps_in,eps_inout_fac,eps_out + real(kind=8) :: ECL,Elj,Equad,Epol,eheadtail,rhead,dGCLOM2,& + sqom1,sqom2,sqom12,c1,c2,c3,pom,Lambf,sparrow,& + Chif,ChiLambf,bat,eagle,top,bot,botsq,Fcav,dtop,dFdR,dFdOM1,& + dFdOM2,w1,w2,w3,dGCLdR,dFdL,dFdOM12,dbot ,& + r1,eps_head,alphapol1,pis,facd2,d2,facd1,d1,erdxj,erdxi,federmaus,& + dPOLdR1,dFGBdOM2,dFGBdR1,dPOLdFGB1,RR1,MomoFac1,hawk,d1i,d1j,& + sig1,sig2,chis12,chis2,ee1,fgb1,a12sq,chis1 + real(kind=8),dimension(3,2)::chead,erhead_tail + real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead + integer troll + eps_out=80.0d0 + epepbase=0.0d0 +! do i=1,nres_molec(1)-1 + do i=ibond_start,ibond_end + if (itype(i,1).eq.ntyp1_molec(1).or.itype(i+1,1).eq.ntyp1_molec(1)) cycle +!C itypi = itype(i,1) + dxi = dc_norm(1,i) + dyi = dc_norm(2,i) + dzi = dc_norm(3,i) +! print *,dxi,(-c(1,i)+c(1,i+1))*vbld_inv(i+1) + dsci_inv = vbld_inv(i+1)/2.0 + 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=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=nres_molec(1)+1,nres_molec(2)+nres_molec(1) + itypj= itype(j,2) + if (itype(j,2).eq.ntyp1_molec(2))cycle + xj=c(1,j+nres) + yj=c(2,j+nres) + zj=c(3,j+nres) + 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 + dxj = dc_norm( 1, nres+j ) + dyj = dc_norm( 2, nres+j ) + dzj = dc_norm( 3, nres+j ) +! d1i = dhead_scbasei(itypi) !this is shift of dipole/charge +! d1j = dhead_scbasej(itypi) !this is shift of dipole/charge + +! Gay-berne var's + sig0ij = sigma_pepbase(itypj ) + chi1 = chi_pepbase(itypj,1 ) + chi2 = chi_pepbase(itypj,2 ) +! chi1=0.0d0 +! chi2=0.0d0 + chi12 = chi1 * chi2 + chip1 = chipp_pepbase(itypj,1 ) + chip2 = chipp_pepbase(itypj,2 ) +! chip1=0.0d0 +! chip2=0.0d0 + chip12 = chip1 * chip2 + chis1 = chis_pepbase(itypj,1) + chis2 = chis_pepbase(itypj,2) + chis12 = chis1 * chis2 + sig1 = sigmap1_pepbase(itypj) + sig2 = sigmap2_pepbase(itypj) +! write (*,*) "sig1 = ", sig1 +! write (*,*) "sig2 = ", sig2 + 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 +! + d1i * dc_norm(k, i+nres) + chead(k,2) = c(k, j+nres) +! + d1j * dc_norm(k, j+nres) +! 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) +! print *,gvdwc_pepbase(k,i) + + END DO + Rhead = dsqrt( & + (Rhead_distance(1)*Rhead_distance(1)) & + + (Rhead_distance(2)*Rhead_distance(2)) & + + (Rhead_distance(3)*Rhead_distance(3))) + +! alpha factors from Fcav/Gcav + b1 = alphasur_pepbase(1,itypj) +! b1=0.0d0 + b2 = alphasur_pepbase(2,itypj) + b3 = alphasur_pepbase(3,itypj) + b4 = alphasur_pepbase(4,itypj) + alf1 = 0.0d0 + alf2 = 0.0d0 + alf12 = 0.0d0 + rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj) +! print *,i,j,rrij + rij = dsqrt(rrij) +!---------------------------- + 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) + 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.0/rij - 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_pepbase(itypj) +! c1 = 0.0d0 + c2 = fac * bb_pepbase(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 + c1 = c1 * eps1 * eps2rt**2 * eps3rt**2 + fac = -expon * (c1 + evdwij) * rij_shift + sigder = fac * sigder +! fac = rij * fac +! Calculate distance derivative + gg(1) = fac + gg(2) = fac + gg(3) = fac + fac = chis1 * sqom1 + chis2 * sqom2 & + - 2.0d0 * chis12 * om1 * om2 * om12 +! we will use pom later in Gcav, so dont mess with it! + pom = 1.0d0 - chis1 * chis2 * sqom12 + Lambf = (1.0d0 - (fac / pom)) + Lambf = dsqrt(Lambf) + sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0) +! write (*,*) "sparrow = ", sparrow + Chif = 1.0d0/rij * sparrow + ChiLambf = Chif * Lambf + eagle = dsqrt(ChiLambf) + bat = ChiLambf ** 11.0d0 + top = b1 * ( eagle + b2 * ChiLambf - b3 ) + bot = 1.0d0 + b4 * (ChiLambf ** 12.0d0) + botsq = bot * bot + Fcav = top / bot +! print *,i,j,Fcav + dtop = b1 * ((Lambf / (2.0d0 * eagle)) + (b2 * Lambf)) + dbot = 12.0d0 * b4 * bat * Lambf + dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow +! dFdR = 0.0d0 +! write (*,*) "dFcav/dR = ", dFdR + dtop = b1 * ((Chif / (2.0d0 * eagle)) + (b2 * Chif)) + dbot = 12.0d0 * b4 * 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) +! dFdL = 0.0d0 + dCAVdOM1 = dFdL * ( dFdOM1 ) + dCAVdOM2 = dFdL * ( dFdOM2 ) + dCAVdOM12 = dFdL * ( dFdOM12 ) + + ertail(1) = xj*rij + ertail(2) = yj*rij + ertail(3) = zj*rij + DO k = 1, 3 +! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i) +! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j) + pom = ertail(k) +!-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres)) + gvdwc_pepbase(k,i) = gvdwc_pepbase(k,i) & + - (( dFdR + gg(k) ) * pom)/2.0 +! print *,gvdwc_pepbase(k,i),i,(( dFdR + gg(k) ) * pom)/2.0 +! +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) & +! +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv +! & - ( dFdR * pom ) + pom = ertail(k) +!-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres)) + gvdwx_pepbase(k,j) = gvdwx_pepbase(k,j) & + + (( dFdR + gg(k) ) * pom) +! +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) & +! +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv +!c! & + ( dFdR * pom ) + + gvdwc_pepbase(k,i+1) = gvdwc_pepbase(k,i+1) & + - (( dFdR + gg(k) ) * ertail(k))/2.0 +! print *,gvdwc_pepbase(k,i+1),i+1,(( dFdR + gg(k) ) * pom)/2.0 + +!c! & - ( dFdR * ertail(k)) + + gvdwc_pepbase(k,j) = gvdwc_pepbase(k,j) & + + (( dFdR + gg(k) ) * ertail(k)) +!c! & + ( dFdR * ertail(k)) + + gg(k) = 0.0d0 +!c! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i) +!c! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j) + END DO + + + w1 = wdipdip_pepbase(1,itypj) + w2 = -wdipdip_pepbase(3,itypj)/2.0 + w3 = wdipdip_pepbase(2,itypj) +! w1=0.0d0 +! w2=0.0d0 +!c!------------------------------------------------------------------- +!c! ECL +! w3=0.0d0 + fac = (om12 - 3.0d0 * om1 * om2) + c1 = (w1 / (Rhead**3.0d0)) * fac + c2 = (w2 / Rhead ** 6.0d0) & + * (4.0d0 + fac * fac -3.0d0 * (sqom1 + sqom2)) + c3= (w3/ Rhead ** 6.0d0) & + * (2.0d0 - 2.0d0*fac*fac +3.0d0*(sqom1 + sqom2)) + + ECL = c1 - c2 + c3 + + c1 = (-3.0d0 * w1 * fac) / (Rhead ** 4.0d0) + c2 = (-6.0d0 * w2) / (Rhead ** 7.0d0) & + * (4.0d0 + fac * fac - 3.0d0 * (sqom1 + sqom2)) + c3= (-6.0d0 * w3) / (Rhead ** 7.0d0) & + * (2.0d0 - 2.0d0*fac*fac +3.0d0*(sqom1 + sqom2)) + + dGCLdR = c1 - c2 + c3 +!c! dECL/dom1 + c1 = (-3.0d0 * w1 * om2 ) / (Rhead**3.0d0) + c2 = (-6.0d0 * w2) / (Rhead**6.0d0) & + * ( om2 * om12 - 3.0d0 * om1 * sqom2 + om1 ) + c3 =(6.0d0*w3/ Rhead ** 6.0d0)*(om1-2.0d0*(fac)*(-om2)) + dGCLdOM1 = c1 - c2 + c3 +!c! dECL/dom2 + c1 = (-3.0d0 * w1 * om1 ) / (Rhead**3.0d0) + c2 = (-6.0d0 * w2) / (Rhead**6.0d0) & + * ( om1 * om12 - 3.0d0 * sqom1 * om2 + om2 ) + c3 =(6.0d0*w3/ Rhead ** 6.0d0)*(om2-2.0d0*(fac)*(-om1)) + + dGCLdOM2 = c1 - c2 + c3 +!c! dECL/dom12 + c1 = w1 / (Rhead ** 3.0d0) + c2 = ( 2.0d0 * w2 * fac ) / Rhead ** 6.0d0 + c3 = (w3/ Rhead ** 6.0d0)*(-4.0d0*fac) + dGCLdOM12 = c1 - c2 + c3 + DO k= 1, 3 + erhead(k) = Rhead_distance(k)/Rhead + END DO + erdxi = scalar( erhead(1), dC_norm(1,i+nres) ) + erdxj = scalar( erhead(1), dC_norm(1,j+nres) ) +! facd1 = d1 * vbld_inv(i+nres) +! facd2 = d2 * vbld_inv(j+nres) + DO k = 1, 3 + +! pom = erhead(k) +!+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres)) +! gvdwx_pepbase(k,i) = gvdwx_scbase(k,i) & +! - dGCLdR * pom + pom = erhead(k) +!+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres)) + gvdwx_pepbase(k,j) = gvdwx_pepbase(k,j) & + + dGCLdR * pom + + gvdwc_pepbase(k,i) = gvdwc_pepbase(k,i) & + - dGCLdR * erhead(k)/2.0d0 +! print *,gvdwc_pepbase(k,i+1),i+1,- dGCLdR * erhead(k)/2.0d0 + gvdwc_pepbase(k,i+1) = gvdwc_pepbase(k,i+1) & + - dGCLdR * erhead(k)/2.0d0 +! print *,gvdwc_pepbase(k,i+1),i+1,- dGCLdR * erhead(k)/2.0d0 + gvdwc_pepbase(k,j) = gvdwc_pepbase(k,j) & + + dGCLdR * erhead(k) + END DO +! print *,i,j,evdwij,Fcav,ECL,"vdw,cav,ecl" + epepbase=epepbase+evdwij+Fcav+ECL + call sc_grad_pepbase + enddo + enddo + END SUBROUTINE epep_sc_base + SUBROUTINE sc_grad_pepbase + use calc_data + + real (kind=8) :: dcosom1(3),dcosom2(3) + 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 +! om12=0.0 +! eom12=0.0 +! print *,eom1,eom2,eom12,om12,i,j,"eom1,2,12",erij(1),erij(2),erij(3) +! if (i.eq.30) print *,gvdwc_pepbase(k,i),- gg(k),& +! (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i)))& +! *dsci_inv*2.0 +! print *,dsci_inv,dscj_inv,dc_norm(2,nres+j),dc_norm(2,nres+i),& +! gg(1),gg(2),"rozne" + 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 +! print *,eom12,eom2,om12,om2 +!eom12*(-dc_norm(k,i)/2.0-om12*dc_norm(k,nres+j)),& +! (eom2*(erij(k)-om2*dc_norm(k,nres+j))) + gvdwx_pepbase(k,j)= gvdwx_pepbase(k,j) + gg(k) & + + (eom12*(dc_norm(k,i)-om12*dc_norm(k,nres+j))& + + eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv + gvdwc_pepbase(k,j)=gvdwc_pepbase(k,j)+gg(k) + END DO + RETURN + END SUBROUTINE sc_grad_pepbase + subroutine eprot_sc_phosphate(escpho) + use calc_data +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +! include 'COMMON.GEO' +! include 'COMMON.VAR' +! include 'COMMON.LOCAL' +! include 'COMMON.CHAIN' +! include 'COMMON.DERIV' +! include 'COMMON.NAMES' +! include 'COMMON.INTERACT' +! include 'COMMON.IOUNITS' +! include 'COMMON.CALC' +! include 'COMMON.CONTROL' +! include 'COMMON.SBRIDGE' + logical :: lprn +!el local variables + integer :: iint,itypi,itypi1,itypj,subchap + real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi + real(kind=8) :: evdw,sig0ij + real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,& + dist_temp, dist_init,aa,bb,ssgradlipi,ssgradlipj, & + sslipi,sslipj,faclip,alpha_sco + integer :: ii + real(kind=8) :: fracinbuf + real (kind=8) :: escpho + real (kind=8),dimension(4):: ener + real(kind=8) :: b1,b2,b3,b4,egb,eps_in,eps_inout_fac,eps_out + real(kind=8) :: ECL,Elj,Equad,Epol,eheadtail,rhead,dGCLOM2,& + sqom1,sqom2,sqom12,c1,c2,pom,Lambf,sparrow,& + Chif,ChiLambf,bat,eagle,top,bot,botsq,Fcav,dtop,dFdR,dFdOM1,& + dFdOM2,w1,w2,dGCLdR,dFdL,dFdOM12,dbot ,& + r1,eps_head,alphapol1,pis,facd2,d2,facd1,d1,erdxj,erdxi,federmaus,& + dPOLdR1,dFGBdOM2,dFGBdR1,dPOLdFGB1,RR1,MomoFac1,hawk,d1i,d1j,& + sig1,sig2,chis12,chis2,ee1,fgb1,a12sq,chis1,Rhead_sq,Qij,dFGBdOM1 + real(kind=8),dimension(3,2)::chead,erhead_tail + real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead + integer troll + eps_out=80.0d0 + escpho=0.0d0 +! do i=1,nres_molec(1) + do i=ibond_start,ibond_end + if (itype(i,1).eq.ntyp1_molec(1)) cycle + itypi = itype(i,1) + dxi = dc_norm(1,nres+i) + dyi = dc_norm(2,nres+i) + dzi = dc_norm(3,nres+i) + dsci_inv = vbld_inv(i+nres) + xi=c(1,nres+i) + yi=c(2,nres+i) + zi=c(3,nres+i) + 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=nres_molec(1)+1,nres_molec(2)+nres_molec(1)-1 + itypj= itype(j,2) + if ((itype(j,2).eq.ntyp1_molec(2)).or.& + (itype(j+1,2).eq.ntyp1_molec(2))) cycle + xj=(c(1,j)+c(1,j+1))/2.0 + yj=(c(2,j)+c(2,j+1))/2.0 + zj=(c(3,j)+c(3,j+1))/2.0 + 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 + dxj = dc_norm( 1,j ) + dyj = dc_norm( 2,j ) + dzj = dc_norm( 3,j ) + dscj_inv = vbld_inv(j+1) + +! Gay-berne var's + sig0ij = sigma_scpho(itypi ) + chi1 = chi_scpho(itypi,1 ) + chi2 = chi_scpho(itypi,2 ) +! chi1=0.0d0 +! chi2=0.0d0 + chi12 = chi1 * chi2 + chip1 = chipp_scpho(itypi,1 ) + chip2 = chipp_scpho(itypi,2 ) +! chip1=0.0d0 +! chip2=0.0d0 + chip12 = chip1 * chip2 + chis1 = chis_scpho(itypi,1) + chis2 = chis_scpho(itypi,2) + chis12 = chis1 * chis2 + sig1 = sigmap1_scpho(itypi) + sig2 = sigmap2_scpho(itypi) +! write (*,*) "sig1 = ", sig1 +! write (*,*) "sig1 = ", sig1 +! write (*,*) "sig2 = ", sig2 +! alpha factors from Fcav/Gcav + alf1 = 0.0d0 + alf2 = 0.0d0 + alf12 = 0.0d0 + a12sq = rborn_scphoi(itypi) * rborn_scphoj(itypi) + + b1 = alphasur_scpho(1,itypi) +! b1=0.0d0 + b2 = alphasur_scpho(2,itypi) + b3 = alphasur_scpho(3,itypi) + b4 = alphasur_scpho(4,itypi) +! used to determine whether we want to do quadrupole calculations +! used by Fgb + eps_in = epsintab_scpho(itypi) + if (eps_in.eq.0.0) eps_in=1.0 + eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out)) +! write (*,*) "eps_inout_fac = ", eps_inout_fac +!------------------------------------------------------------------- +! tail location and distance calculations + d1i = dhead_scphoi(itypi) !this is shift of dipole/charge + d1j = 0.0 + 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+nres) + d1i * dc_norm(k, i+nres) + chead(k,2) = (c(k, j) + c(k, j+1))/2.0 +! 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))) + Rhead_sq=Rhead**2.0 +!------------------------------------------------------------------- +! 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 + dGCLdR=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+1)/2.0 +!dhead_scbasej(itypi,itypj) +! 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 = 1.0/rij - 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_scpho(itypi) +! c1 = 0.0d0 + c2 = fac * bb_scpho(itypi) +! 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 + c1 = c1 * eps1 * eps2rt**2 * eps3rt**2 + fac = -expon * (c1 + evdwij) * rij_shift + sigder = fac * sigder +! fac = rij * fac +! Calculate distance derivative + gg(1) = fac + gg(2) = fac + gg(3) = fac + fac = chis1 * sqom1 + chis2 * sqom2 & + - 2.0d0 * chis12 * om1 * om2 * om12 +! we will use pom later in Gcav, so dont mess with it! + pom = 1.0d0 - chis1 * chis2 * sqom12 + Lambf = (1.0d0 - (fac / pom)) + Lambf = dsqrt(Lambf) + sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0) +! write (*,*) "sparrow = ", sparrow + Chif = 1.0d0/rij * sparrow + ChiLambf = Chif * Lambf + eagle = dsqrt(ChiLambf) + bat = ChiLambf ** 11.0d0 + top = b1 * ( eagle + b2 * ChiLambf - b3 ) + bot = 1.0d0 + b4 * (ChiLambf ** 12.0d0) + botsq = bot * bot + Fcav = top / bot + dtop = b1 * ((Lambf / (2.0d0 * eagle)) + (b2 * Lambf)) + dbot = 12.0d0 * b4 * bat * Lambf + dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow +! dFdR = 0.0d0 +! write (*,*) "dFcav/dR = ", dFdR + dtop = b1 * ((Chif / (2.0d0 * eagle)) + (b2 * Chif)) + dbot = 12.0d0 * b4 * 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) +! dFdL = 0.0d0 + dCAVdOM1 = dFdL * ( dFdOM1 ) + dCAVdOM2 = dFdL * ( dFdOM2 ) + dCAVdOM12 = dFdL * ( dFdOM12 ) + + ertail(1) = xj*rij + ertail(2) = yj*rij + ertail(3) = zj*rij + DO k = 1, 3 +! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i) +! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j) +! if (i.eq.3) print *,'decl0',gvdwx_scpho(k,i),i + + pom = ertail(k) +! print *,pom,gg(k),dFdR +!-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres)) + gvdwx_scpho(k,i) = gvdwx_scpho(k,i) & + - (( dFdR + gg(k) ) * pom) +! +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) & +! +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv +! & - ( dFdR * pom ) +! pom = ertail(k) +!-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres)) +! gvdwx_scpho(k,j) = gvdwx_scpho(k,j) & +! + (( dFdR + gg(k) ) * pom) +! +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) & +! +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv +!c! & + ( dFdR * pom ) + + gvdwc_scpho(k,i) = gvdwc_scpho(k,i) & + - (( dFdR + gg(k) ) * ertail(k)) +!c! & - ( dFdR * ertail(k)) + + gvdwc_scpho(k,j) = gvdwc_scpho(k,j) & + + (( dFdR + gg(k) ) * ertail(k))/2.0 + + gvdwc_scpho(k,j+1) = gvdwc_scpho(k,j+1) & + + (( dFdR + gg(k) ) * ertail(k))/2.0 + +!c! & + ( dFdR * ertail(k)) + + gg(k) = 0.0d0 + ENDDO +!c! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i) +!c! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j) +! alphapol1 = alphapol_scpho(itypi) + if (wqq_scpho(itypi).ne.0.0) then + Qij=wqq_scpho(itypi)/eps_in + alpha_sco=1.d0/alphi_scpho(itypi) +! Qij=0.0 + Ecl = (332.0d0 * Qij*dexp(-Rhead*alpha_sco)) / Rhead +!c! derivative of Ecl is Gcl... + dGCLdR = (-332.0d0 * Qij*dexp(-Rhead*alpha_sco)* & + (Rhead*alpha_sco+1) ) / Rhead_sq + if (energy_dec) write(iout,*) "ECL",ECL,Rhead,1.0/rij + else if (wqdip_scpho(2,itypi).gt.0.0d0) then + w1 = wqdip_scpho(1,itypi) + w2 = wqdip_scpho(2,itypi) +! w1=0.0d0 +! w2=0.0d0 +! pis = sig0head_scbase(itypi,itypj) +! eps_head = epshead_scbase(itypi,itypj) +!c!------------------------------------------------------------------- + +!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 * om1 + hawk = w2 * (1.0d0 - sqom2) + Ecl = sparrow / Rhead**2.0d0 & + - hawk / Rhead**4.0d0 +!c!------------------------------------------------------------------- + if (energy_dec) write(iout,*) "ECLdipdip",ECL,Rhead,& + 1.0/rij,sparrow + +!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) / (Rhead**2.0d0) +!c! dF/dom2 + dGCLdOM2 = (2.0d0 * w2 * om2) / (Rhead ** 4.0d0) + endif + +!c-------------------------------------------------------------------- +!c Polarization energy +!c Epol + R1 = 0.0d0 + DO k = 1, 3 +!c! Calculate head-to-tail distances tail is center of side-chain + R1=R1+((c(k,j)+c(k,j+1))/2.0-chead(k,1))**2 + END DO +!c! Pitagoras + R1 = dsqrt(R1) + + alphapol1 = alphapol_scpho(itypi) +! alphapol1=0.0 + MomoFac1 = (1.0d0 - chi2 * sqom1) + RR1 = R1 * R1 / MomoFac1 + ee1 = exp(-( RR1 / (4.0d0 * a12sq) )) +! print *,"ee1",ee1,a12sq,alphapol1,eps_inout_fac + fgb1 = sqrt( RR1 + a12sq * ee1) +! eps_inout_fac=0.0d0 + epol = 332.0d0 * eps_inout_fac * (( alphapol1 / fgb1 )**4.0d0) +! derivative of Epol is Gpol... + dPOLdFGB1 = -(1328.0d0 * eps_inout_fac * alphapol1 ** 4.0d0) & + / (fgb1 ** 5.0d0) + dFGBdR1 = ( (R1 / MomoFac1) & + * ( 2.0d0 - (0.5d0 * ee1) ) ) & + / ( 2.0d0 * fgb1 ) + dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1)) & + * (2.0d0 - 0.5d0 * ee1) ) & + / (2.0d0 * fgb1) + dPOLdR1 = dPOLdFGB1 * dFGBdR1 +! dPOLdR1 = 0.0d0 +! dPOLdOM1 = 0.0d0 + dFGBdOM1 = (((R1 * R1 * chi2 * om1) / (MomoFac1 * MomoFac1)) & + * (2.0d0 - 0.5d0 * ee1) ) & + / (2.0d0 * fgb1) + + dPOLdOM1 = dPOLdFGB1 * dFGBdOM1 + dPOLdOM2 = 0.0 + DO k = 1, 3 + erhead(k) = Rhead_distance(k)/Rhead + erhead_tail(k,1) = (((c(k,j)+c(k,j+1))/2.0-chead(k,1))/R1) + END DO + + erdxi = scalar( erhead(1), dC_norm(1,i+nres) ) + erdxj = scalar( erhead(1), dC_norm(1,j) ) + bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) ) +! bat=0.0d0 + federmaus = scalar(erhead_tail(1,1),dC_norm(1,j)) + facd1 = d1i * vbld_inv(i+nres) + facd2 = d1j * vbld_inv(j) +! facd4 = dtail(2,itypi,itypj) * vbld_inv(j+nres) + + DO k = 1, 3 + hawk = (erhead_tail(k,1) + & + facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres))) +! facd1=0.0d0 +! facd2=0.0d0 +! if (i.eq.3) print *,'decl1',dGCLdR,dPOLdR1,gvdwc_scpho(k,i),i,& +! pom,(erhead_tail(k,1)) + +! print *,'decl',dGCLdR,dPOLdR1,gvdwc_scpho(k,i) + pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres)) + gvdwx_scpho(k,i) = gvdwx_scpho(k,i) & + - dGCLdR * pom & + - dPOLdR1 * (erhead_tail(k,1)) +! & - dGLJdR * pom + + pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j)) +! gvdwx_scpho(k,j) = gvdwx_scpho(k,j) & +! + dGCLdR * pom & +! + dPOLdR1 * (erhead_tail(k,1)) +! & + dGLJdR * pom + + + gvdwc_scpho(k,i) = gvdwc_scpho(k,i) & + - dGCLdR * erhead(k) & + - dPOLdR1 * erhead_tail(k,1) +! & - dGLJdR * erhead(k) + + gvdwc_scpho(k,j) = gvdwc_scpho(k,j) & + + (dGCLdR * erhead(k) & + + dPOLdR1 * erhead_tail(k,1))/2.0 + gvdwc_scpho(k,j+1) = gvdwc_scpho(k,j+1) & + + (dGCLdR * erhead(k) & + + dPOLdR1 * erhead_tail(k,1))/2.0 + +! & + dGLJdR * erhead(k) +! if (i.eq.3) print *,'decl2',dGCLdR,dPOLdR1,gvdwc_scpho(k,i),i + + END DO +! if (i.eq.3) print *,i,j,evdwij,epol,Fcav,ECL + if (energy_dec) write (iout,'(a22,2i5,4f8.3,f16.3)'), & + "escpho:evdw,pol,cav,CL",i,j,evdwij,epol,Fcav,ECL,escpho + escpho=escpho+evdwij+epol+Fcav+ECL + call sc_grad_scpho + enddo + + enddo + + return + end subroutine eprot_sc_phosphate + SUBROUTINE sc_grad_scpho + use calc_data + + real (kind=8) :: dcosom1(3),dcosom2(3) + 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 +! om12=0.0 +! eom12=0.0 +! print *,eom1,eom2,eom12,om12,i,j,"eom1,2,12",erij(1),erij(2),erij(3) +! if (i.eq.30) print *,gvdwc_scpho(k,i),- gg(k),& +! (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i)))& +! *dsci_inv*2.0 +! print *,dsci_inv,dscj_inv,dc_norm(2,nres+j),dc_norm(2,nres+i),& +! gg(1),gg(2),"rozne" + DO k = 1, 3 + dcosom1(k) = rij * (dc_norm(k,nres+i) - om1 * erij(k)) + dcosom2(k) = rij * (dc_norm(k,j) - om2 * erij(k)) + gg(k) = gg(k) + eom1 * dcosom1(k) + eom2 * dcosom2(k) + gvdwc_scpho(k,j)= gvdwc_scpho(k,j) +0.5*( gg(k)) & + + (-eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,j)))& + *dscj_inv*2.0 & + - (eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv*2.0 + gvdwc_scpho(k,j+1)= gvdwc_scpho(k,j+1) +0.5*( gg(k)) & + - (-eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,j))) & + *dscj_inv*2.0 & + + (eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv*2.0 + gvdwx_scpho(k,i)= gvdwx_scpho(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 + +! print *,eom12,eom2,om12,om2 +!eom12*(-dc_norm(k,i)/2.0-om12*dc_norm(k,nres+j)),& +! (eom2*(erij(k)-om2*dc_norm(k,nres+j))) +! gvdwx_scpho(k,j)= gvdwx_scpho(k,j) + gg(k) & +! + (eom12*(dc_norm(k,i)-om12*dc_norm(k,nres+j))& +! + eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv + gvdwc_scpho(k,i)=gvdwc_scpho(k,i)-gg(k) + END DO + RETURN + END SUBROUTINE sc_grad_scpho + subroutine eprot_pep_phosphate(epeppho) + use calc_data +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +! include 'COMMON.GEO' +! include 'COMMON.VAR' +! include 'COMMON.LOCAL' +! include 'COMMON.CHAIN' +! include 'COMMON.DERIV' +! include 'COMMON.NAMES' +! include 'COMMON.INTERACT' +! include 'COMMON.IOUNITS' +! include 'COMMON.CALC' +! include 'COMMON.CONTROL' +! include 'COMMON.SBRIDGE' + logical :: lprn +!el local variables + integer :: iint,itypi,itypi1,itypj,subchap + real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi + real(kind=8) :: evdw,sig0ij + real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,& + dist_temp, dist_init,aa,bb,ssgradlipi,ssgradlipj, & + sslipi,sslipj,faclip + integer :: ii + real(kind=8) :: fracinbuf + real (kind=8) :: epeppho + real (kind=8),dimension(4):: ener + real(kind=8) :: b1,b2,b3,b4,egb,eps_in,eps_inout_fac,eps_out + real(kind=8) :: ECL,Elj,Equad,Epol,eheadtail,rhead,dGCLOM2,& + sqom1,sqom2,sqom12,c1,c2,pom,Lambf,sparrow,& + Chif,ChiLambf,bat,eagle,top,bot,botsq,Fcav,dtop,dFdR,dFdOM1,& + dFdOM2,w1,w2,dGCLdR,dFdL,dFdOM12,dbot ,& + r1,eps_head,alphapol1,pis,facd2,d2,facd1,d1,erdxj,erdxi,federmaus,& + dPOLdR1,dFGBdOM2,dFGBdR1,dPOLdFGB1,RR1,MomoFac1,hawk,d1i,d1j,& + sig1,sig2,chis12,chis2,ee1,fgb1,a12sq,chis1,Rhead_sq,Qij,dFGBdOM1 + real(kind=8),dimension(3,2)::chead,erhead_tail + real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead + integer troll + real (kind=8) :: dcosom1(3),dcosom2(3) + epeppho=0.0d0 +! do i=1,nres_molec(1) + do i=ibond_start,ibond_end + if (itype(i,1).eq.ntyp1_molec(1)) cycle + itypi = itype(i,1) + dsci_inv = vbld_inv(i+1)/2.0 + dxi = dc_norm(1,i) + dyi = dc_norm(2,i) + dzi = dc_norm(3,i) + 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=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=nres_molec(1)+1,nres_molec(2)+nres_molec(1)-1 + itypj= itype(j,2) + if ((itype(j,2).eq.ntyp1_molec(2)).or.& + (itype(j+1,2).eq.ntyp1_molec(2))) cycle + xj=(c(1,j)+c(1,j+1))/2.0 + yj=(c(2,j)+c(2,j+1))/2.0 + zj=(c(3,j)+c(3,j+1))/2.0 + 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 + rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj) + rij = dsqrt(rrij) + dxj = dc_norm( 1,j ) + dyj = dc_norm( 2,j ) + dzj = dc_norm( 3,j ) + dscj_inv = vbld_inv(j+1)/2.0 +! Gay-berne var's + sig0ij = sigma_peppho + chi1=0.0d0 + chi2=0.0d0 + chi12 = chi1 * chi2 + chip1=0.0d0 + chip2=0.0d0 + chip12 = chip1 * chip2 + chis1 = 0.0d0 + chis2 = 0.0d0 + chis12 = chis1 * chis2 + sig1 = sigmap1_peppho + sig2 = sigmap2_peppho +! write (*,*) "sig1 = ", sig1 +! write (*,*) "sig1 = ", sig1 +! write (*,*) "sig2 = ", sig2 +! alpha factors from Fcav/Gcav + alf1 = 0.0d0 + alf2 = 0.0d0 + alf12 = 0.0d0 + b1 = alphasur_peppho(1) +! b1=0.0d0 + b2 = alphasur_peppho(2) + b3 = alphasur_peppho(3) + b4 = alphasur_peppho(4) + CALL sc_angular + sqom1=om1*om1 + evdwij = 0.0d0 + ECL = 0.0d0 + Elj = 0.0d0 + Equad = 0.0d0 + Epol = 0.0d0 + Fcav=0.0d0 + eheadtail = 0.0d0 + dGCLdR=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 + rij_shift = rij + fac = rij_shift**expon + c1 = fac * fac * aa_peppho +! c1 = 0.0d0 + c2 = fac * bb_peppho +! c2 = 0.0d0 + evdwij = c1 + c2 +! Now cavity.................... + eagle = dsqrt(1.0/rij_shift) + top = b1 * ( eagle + b2 * 1.0/rij_shift - b3 ) + bot = 1.0d0 + b4 * (1.0/rij_shift ** 12.0d0) + botsq = bot * bot + Fcav = top / bot + dtop = b1 * ((1.0/ (2.0d0 * eagle)) + (b2)) + dbot = 12.0d0 * b4 * (1.0/rij_shift) ** 11.0d0 + dFdR = ((dtop * bot - top * dbot) / botsq) + w1 = wqdip_peppho(1) + w2 = wqdip_peppho(2) +! w1=0.0d0 +! w2=0.0d0 +! pis = sig0head_scbase(itypi,itypj) +! eps_head = epshead_scbase(itypi,itypj) +!c!------------------------------------------------------------------- + +!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 * om1 + hawk = w2 * (1.0d0 - sqom1) + Ecl = sparrow * rij_shift**2.0d0 & + - hawk * rij_shift**4.0d0 +!c!------------------------------------------------------------------- +!c! derivative of ecl is Gcl +!c! dF/dr part +! rij_shift=5.0 + dGCLdR = - 2.0d0 * sparrow * rij_shift**3.0d0 & + + 4.0d0 * hawk * rij_shift**5.0d0 +!c! dF/dom1 + dGCLdOM1 = (w1) * (rij_shift**2.0d0) +!c! dF/dom2 + dGCLdOM2 = (2.0d0 * w2 * om1) * (rij_shift ** 4.0d0) + eom1 = dGCLdOM1+dGCLdOM2 + eom2 = 0.0 + + fac = -expon * (c1 + evdwij) * rij_shift+dFdR+dGCLdR +! fac=0.0 + gg(1) = fac*xj*rij + gg(2) = fac*yj*rij + gg(3) = fac*zj*rij + do k=1,3 + gvdwc_peppho(k,j) = gvdwc_peppho(k,j) +gg(k)/2.0 + gvdwc_peppho(k,j+1) = gvdwc_peppho(k,j+1) +gg(k)/2.0 + gvdwc_peppho(k,i) = gvdwc_peppho(k,i) -gg(k)/2.0 + gvdwc_peppho(k,i+1) = gvdwc_peppho(k,i+1) -gg(k)/2.0 + gg(k)=0.0 + enddo + + DO k = 1, 3 + dcosom1(k) = rij* (dc_norm(k,i) - om1 * erij(k)) + dcosom2(k) = rij* (dc_norm(k,j) - om2 * erij(k)) + gg(k) = gg(k) + eom1 * dcosom1(k)! + eom2 * dcosom2(k) + gvdwc_peppho(k,j)= gvdwc_peppho(k,j) +0.5*( gg(k)) !& +! - (eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv*2.0 + gvdwc_peppho(k,j+1)= gvdwc_peppho(k,j+1) +0.5*( gg(k)) !& +! + (eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv*2.0 + gvdwc_peppho(k,i)= gvdwc_peppho(k,i) -0.5*( gg(k)) & + - (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0 + gvdwc_peppho(k,i+1)= gvdwc_peppho(k,i+1) - 0.5*( gg(k)) & + + (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0 + enddo + epeppho=epeppho+evdwij+Fcav+ECL +! print *,i,j,evdwij,Fcav,ECL,rij_shift + enddo + enddo + end subroutine eprot_pep_phosphate end module energy