X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Funres%2Fenergy.f90;h=ac1202871443616851aa948f8c764d29a98e4a6f;hb=dc956ab7ea5a327e869c9671337d951a55c2535b;hp=c92f4bc948c4afc8c34d3a6f15dfdc07fdfe0075;hpb=254aac5de5a2d56e17a059d77c5f32f7b77a41fa;p=unres4.git diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index c92f4bc..ac12028 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -130,10 +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, & 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,& @@ -241,6 +247,8 @@ 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 @@ -294,8 +302,11 @@ weights_(36)=wtor_d_nucl weights_(37)=wcorr_nucl weights_(38)=wcorr3_nucl - weights(41)=wcatcat - weights(42)=wcatprot + weights_(41)=wcatcat + weights_(42)=wcatprot + weights_(46)=wscbase + weights_(47)=wscpho + weights_(48)=wpeppho ! wcatcat= weights(41) ! wcatprot=weights(42) @@ -340,6 +351,9 @@ 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 @@ -601,6 +615,10 @@ call multibody_hb_nucl(ecorr_nucl,ecorr3_nucl,n_corr,n_corr1) call ecatcat(ecationcation) call ecat_prot(ecation_prot) + call eprot_sc_base(escbase) + call epep_sc_base(epepbase) + call eprot_sc_phosphate(escpho) + call eprot_pep_phosphate(epeppho) ! call ecatcat(ecationcation) ! print *,"after ebend", ebe_nucl #ifdef TIMING @@ -668,6 +686,10 @@ ! 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" @@ -711,6 +733,7 @@ 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 @@ -789,6 +812,10 @@ 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 @@ -805,7 +832,8 @@ +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& - +wcatprot*ecation_prot+wcatcat*ecationcation + +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 & @@ -818,7 +846,8 @@ +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& - +wcatprot*ecation_prot+wcatcat*ecationcation + +wcatprot*ecation_prot+wcatcat*ecationcation+wscbase*escbase& + +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho #endif energia(0)=etot ! detecting NaNQ @@ -947,6 +976,7 @@ 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) @@ -995,7 +1025,10 @@ 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,& @@ -1010,6 +1043,7 @@ evdwsb,wvdwsb,eelsb,welsb,esbloc,wsbloc,etors_nucl,wtor_nucl,& 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,& etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & @@ -1053,6 +1087,10 @@ '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,& @@ -1068,6 +1106,7 @@ evdwsb,wvdwsb,eelsb,welsb,esbloc,wsbloc,etors_nucl,wtor_nucl& 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,& etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & @@ -1110,6 +1149,10 @@ '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 @@ -10476,7 +10519,14 @@ wcorr_nucl*gradcorr_nucl(j,i)& +wcorr3_nucl*gradcorr3_nucl(j,i)+& wcatprot* gradpepcat(j,i)+ & - wcatcat*gradcatcat(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 @@ -10507,7 +10557,13 @@ wcorr_nucl*gradcorr_nucl(j,i) & +wcorr3_nucl*gradcorr3_nucl(j,i) +& wcatprot* gradpepcat(j,i)+ & - wcatcat*gradcatcat(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 @@ -10766,7 +10822,12 @@ +wcorr_nucl*gradxorr_nucl(j,i)& +wcorr3_nucl*gradxorr3_nucl(j,i) & +wsbloc*gsblocx(j,i) & - +wcatprot* gradpepcatx(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 @@ -16373,6 +16434,13 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -19783,6 +19851,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)) @@ -22293,4 +22369,1692 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !---------------------------------------------------------------------------- !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- + 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 + 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 + 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).gt.0.0) then + Qij=wqq_scpho(itypi)/eps_in +! Qij=0.0 + Ecl = (332.0d0 * Qij) / Rhead +!c! derivative of Ecl is Gcl... + dGCLdR = (-332.0d0 * Qij ) / Rhead_sq + 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!------------------------------------------------------------------- +!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 + 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 + 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