From e42cb389c07a8bdb6de95554720f33e09c701cce Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Wed, 29 Nov 2017 10:36:42 +0100 Subject: [PATCH] working peptide-phosphate --- source/unres/MREMD.f90 | 10 +- source/unres/control.F90 | 2 + source/unres/data/energy_data.f90 | 7 +- source/unres/data/io_units.f90 | 4 +- source/unres/data/names.f90 | 8 +- source/unres/energy.f90 | 270 +++++++++++++++++++++++++++++++++++-- source/unres/io.f90 | 2 + source/unres/io_config.f90 | 21 +++ 8 files changed, 300 insertions(+), 24 deletions(-) diff --git a/source/unres/MREMD.f90 b/source/unres/MREMD.f90 index 350ab52..97a91e5 100644 --- a/source/unres/MREMD.f90 +++ b/source/unres/MREMD.f90 @@ -84,7 +84,7 @@ real(kind=8) :: energia(0:n_ene) real(kind=8) :: remd_t_bath(Nprocs) !(maxprocs) integer :: iremd_iset(Nprocs) !(maxprocs) - integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200) + integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,max0(Nprocs/200,1),max0(Nprocs/200,1)) ! (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200) real(kind=8) :: remd_ene(0:n_ene+4,Nprocs) !(0:n_ene+4,maxprocs) integer :: iremd_acc(Nprocs),iremd_tot(Nprocs) !(maxprocs) @@ -1341,7 +1341,9 @@ real(kind=4) :: r_d(3,0:2*nres) real(kind=4) :: t5_restart1(5) integer :: iret,itmp - integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200) +! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200) + integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,max0(Nprocs/200,1),max0(Nprocs/200,1)) + !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200) !el common /przechowalnia/ d_restart1,d_restart2 integer :: i,j,il,il1,ierr,ixdrf @@ -1755,7 +1757,9 @@ ! include 'COMMON.INTERACT' !el real(kind=4) :: d_restart1(3,2*nres*maxprocs) real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5) - integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200) +! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200) + integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,max0(Nprocs/200,1),max0(Nprocs/200,1)) + !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200) !el common /przechowalnia/ d_restart1 integer :: i,j,il,il1,ierr,itmp,iret,ixdrf diff --git a/source/unres/control.F90 b/source/unres/control.F90 index 55d9f9c..0444b44 100644 --- a/source/unres/control.F90 +++ b/source/unres/control.F90 @@ -240,6 +240,8 @@ isidep_scbase=141 isidep_pepbase=142 isidep_scpho=143 + isidep_peppho=144 + iliptranpar=60 itube=61 ! IONS diff --git a/source/unres/data/energy_data.f90 b/source/unres/data/energy_data.f90 index bb00815..6083686 100644 --- a/source/unres/data/energy_data.f90 +++ b/source/unres/data/energy_data.f90 @@ -67,7 +67,7 @@ wturn6,wvdwpp,wliptran,wshield,lipscale,wtube, & wbond_nucl,wang_nucl,wcorr_nucl,wcorr3_nucl,welpp,wtor_nucl,& wtor_d_nucl,welsb,wsbloc,wvdwsb,welpsb,wvdwpp_nucl,wvdwpsb,wcatprot,& - wcatcat,wscbase,wpepbase,wscpho + wcatcat,wscbase,wpepbase,wscpho,wpeppho #ifdef CLUSTER real(kind=8) :: scalscp #endif @@ -365,5 +365,8 @@ real(kind=8),dimension(:,:), allocatable :: alphasur_scpho, & chi_scpho,chipp_scpho,chis_scpho, & wqdip_scpho - + real(kind=8) ,dimension(4) :: alphasur_peppho + real(kind=8) ,dimension(2) :: wqdip_peppho + real(kind=8) :: eps_peppho,sigma_peppho,sigmap1_peppho,sigmap2_peppho, & + aa_peppho,bb_peppho end module energy_data diff --git a/source/unres/data/io_units.f90 b/source/unres/data/io_units.f90 index 762eda0..909e6fd 100644 --- a/source/unres/data/io_units.f90 +++ b/source/unres/data/io_units.f90 @@ -17,7 +17,7 @@ ithep_pdb,irotam_pdb,iliptranpar,itube, & ibond_nucl,ithep_nucl,irotam_nucl,itorp_nucl,itordp_nucl, & isidep_nucl,iscpp_nucl,isidep_scbase,isidep_pepbase,isidep_scpho,& - iion + isidep_peppho,iion #ifdef WHAM_RUN ! el wham iounits integer :: isidep1,ihist,iweight,izsc,idistr @@ -49,7 +49,7 @@ thetname_pdb,rotname_pdb,liptranname,tubename, & bondname_nucl,thetname_nucl,rotname_nucl,torname_nucl, & tordname_nucl,sidename_nucl,scpname_nucl, & - sidename_scbase,pepname_pepbase,pepname_scpho, & + sidename_scbase,pepname_pepbase,pepname_scpho,pepname_peppho, & ionname !----------------------------------------------------------------------- ! INP - main input file diff --git a/source/unres/data/names.f90 b/source/unres/data/names.f90 index 7b26e0f..dc8b2e1 100644 --- a/source/unres/data/names.f90 +++ b/source/unres/data/names.f90 @@ -59,7 +59,7 @@ !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! Number of energy components - integer,parameter :: n_ene=48 + integer,parameter :: n_ene=49 integer :: n_ene2=2*n_ene !----------------------------------------------------------------------------- ! common.names @@ -88,7 +88,7 @@ "EESSB ","ESTR ","EBE ","ESBLOC ","ETORS ",& "ETORSD ","ECORR ","ECORR3 ","NULL ","NULL ",& "ECATPROT ","ECATCAT ","NULL ","NULL ","NULL ",& - "ESCBASE ","EPEPBASE ","ESCPHO "/) + "ESCBASE ","EPEPBASE ","ESCPHO ","EPEPPHO "/) character(len=10),dimension(n_ene) :: wname = & (/"WSC ","WSCP ","WELEC" ,"WCORR ","WCORR5 ","WCORR6 ","WEL_LOC ",& @@ -99,13 +99,13 @@ "WELSB ","WBOND ","WANG ","WSBLOC ","WTOR ",& "WTORD ","WCORR ","WCORR3 ","WNULL ","WNULL ",& "WCATPROT ","WCATCAT ","WNULL ","WNULL ","WNULL ",& - "WSCBASE ","WPEPBASE ","WSCPHO "/) + "WSCBASE ","WPEPBASE ","WSCPHO ","WPEPPH O "/) integer :: nprint_ene = 21 integer,dimension(n_ene) :: print_order = & (/1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,16,15,17,20,21,22,23,24,25,& 26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,& - 48/) + 48,49/) character(len=1), dimension(2) :: sugartyp = (/'D',' '/) !#endif diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index 4efd69f..a3b45d0 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -132,7 +132,8 @@ gvdwpp_nucl !-----------------------------NUCLEIC-PROTEIN GRADIENT real(kind=8),dimension(:,:),allocatable :: gvdwx_scbase,gvdwc_scbase,& - gvdwx_pepbase,gvdwc_pepbase,gvdwx_scpho,gvdwc_scpho + gvdwx_pepbase,gvdwc_pepbase,gvdwx_scpho,gvdwc_scpho,& + gvdwc_peppho !------------------------------IONS GRADIENT real(kind=8),dimension(:,:),allocatable :: gradcatcat, & gradpepcat,gradpepcatx @@ -247,7 +248,7 @@ ! energies for ions real(kind=8) :: ecation_prot,ecationcation ! energies for protein nucleic acid interaction - real(kind=8) :: escbase,epepbase,escpho + real(kind=8) :: escbase,epepbase,escpho,epeppho #ifdef MPI real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw @@ -305,6 +306,7 @@ weights_(42)=wcatprot weights_(46)=wscbase weights_(47)=wscpho + weights_(48)=wpeppho ! wcatcat= weights(41) ! wcatprot=weights(42) @@ -351,6 +353,7 @@ 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 @@ -615,6 +618,7 @@ 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 @@ -685,6 +689,7 @@ 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" @@ -728,7 +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 + real(kind=8) :: escbase,epepbase,escpho,epeppho integer :: i #ifdef MPI integer :: ierr @@ -810,6 +815,7 @@ escbase=energia(46) epepbase=energia(47) escpho=energia(48) + epeppho=energia(49) ! energia(41)=ecation_prot ! energia(42)=ecationcation @@ -827,7 +833,7 @@ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl& +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl& +wcatprot*ecation_prot+wcatcat*ecationcation+wscbase*escbase& - +wpepbase*epepbase+wscpho*escpho + +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc & @@ -841,7 +847,7 @@ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl& +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl& +wcatprot*ecation_prot+wcatcat*ecationcation+wscbase*escbase& - +wpepbase*epepbase+wscpho*escpho + +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho #endif energia(0)=etot ! detecting NaNQ @@ -970,7 +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 + real(kind=8) :: escbase,epepbase,escpho,epeppho etot=energia(0) evdw=energia(1) @@ -1022,6 +1028,7 @@ 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,& @@ -1036,7 +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,& + escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,& etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & @@ -1083,6 +1090,7 @@ '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,& @@ -1098,7 +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,& + escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,& etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & @@ -1144,6 +1152,7 @@ '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 @@ -10513,7 +10522,9 @@ wcatcat*gradcatcat(j,i)+ & wscbase*gvdwc_scbase(j,i)+ & wpepbase*gvdwc_pepbase(j,i)+& - wscpho*gvdwc_scpho(j,i) + wscpho*gvdwc_scpho(j,i)+ & + wpeppho*gvdwc_peppho(j,i) + @@ -10549,7 +10560,8 @@ wcatcat*gradcatcat(j,i)+ & wscbase*gvdwc_scbase(j,i) & wpepbase*gvdwc_pepbase(j,i)+& - wscpho*gvdwc_scpho(j,i) + wscpho*gvdwc_scpho(j,i)+& + wpeppho*gvdwc_peppho(j,i) enddo @@ -10814,7 +10826,7 @@ +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) +! if (i.eq.3) print *,"tu?", wscpho,gvdwx_scpho(j,i) enddo enddo @@ -16428,6 +16440,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -19844,6 +19857,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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) @@ -22671,7 +22685,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !Now dipole-dipole if (wdipdip_scbase(2,itypi,itypj).gt.0.0d0) then w1 = wdipdip_scbase(1,itypi,itypj) - w2 = wdipdip_scbase(2,itypi,itypj) + w2 = wdipdip_scbase(3,itypi,itypj)*2.0 !c!------------------------------------------------------------------- !c! ECL fac = (om12 - 3.0d0 * om1 * om2) @@ -23171,7 +23185,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' w1 = wdipdip_pepbase(1,itypj) - w2 = wdipdip_pepbase(2,itypj) + w2 = wdipdip_pepbase(3,itypj)*2.0 ! w1=0.0d0 ! w2=0.0d0 !c!------------------------------------------------------------------- @@ -23391,6 +23405,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 ) @@ -23787,5 +23803,233 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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) + 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 diff --git a/source/unres/io.f90 b/source/unres/io.f90 index 5ecb9a7..0458733 100644 --- a/source/unres/io.f90 +++ b/source/unres/io.f90 @@ -803,6 +803,8 @@ call reada(weightcard,'WCATPROT',wcatprot,0.0d0) call reada(weightcard,'WPEPBASE',wpepbase,1.0d0) call reada(weightcard,'WSCPHO',wscpho,1.0d0) + call reada(weightcard,'WPEPPHO',wpeppho,1.0d0) + ! 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WCORRH',wcorr,1.0D0) if (wcorr4.gt.0.0d0) wcorr=wcorr4 diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90 index 2d4a625..91953ab 100644 --- a/source/unres/io_config.f90 +++ b/source/unres/io_config.f90 @@ -2550,6 +2550,25 @@ enddo + read(isidep_peppho,*) & + eps_peppho,sigma_peppho + read(isidep_peppho,*) & + (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho + read(isidep_peppho,*) & + (wqdip_peppho(k),k=1,2) + + epsij=eps_peppho + rrij=sigma_peppho +! r0(i,j)=rrij +! r0(j,i)=rrij + rrij=rrij**expon +! epsij=eps(i,j) + sigeps=dsign(1.0D0,epsij) + epsij=dabs(epsij) + aa_peppho=epsij*rrij*rrij + bb_peppho=-sigeps*epsij*rrij + + allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2) bad(:,:)=0.0D0 @@ -4554,6 +4573,8 @@ open (isidep_pepbase,file=pepname_pepbase,status='old',action='read') call getenv_loc('SCPAR_PHOSPH',pepname_scpho) open (isidep_scpho,file=pepname_scpho,status='old',action='read') + call getenv_loc('PEPPAR_PHOSPH',pepname_peppho) + open (isidep_peppho,file=pepname_peppho,status='old',action='read') call getenv_loc('LIPTRANPAR',liptranname) -- 1.7.9.5