X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Funres%2Fenergy.F90;h=5c9aacf99464db97144f92f75ab04ce394a5046f;hb=f8ee39d07f69eabbbf4945b63c93308c85fe5a96;hp=95ebd9d567e246cf10e00760ac9c7cec63699149;hpb=df2469d9ac903d93889867f4e50e9bf6c428c1c6;p=unres4.git diff --git a/source/unres/energy.F90 b/source/unres/energy.F90 index 95ebd9d..5c9aacf 100644 --- a/source/unres/energy.F90 +++ b/source/unres/energy.F90 @@ -1,4 +1,4 @@ - module energy + module energy !----------------------------------------------------------------------------- use io_units use names @@ -74,7 +74,7 @@ ! amino-acid residue. ! common /precomp1/ real(kind=8),dimension(:,:),allocatable :: mu,muder,Ub2,Ub2der,& - Ctobr,Ctobrder,Dtobr2,Dtobr2der !(2,maxres) + Ctobr,Ctobrder,Dtobr2,Dtobr2der,gUb2 !(2,maxres) real(kind=8),dimension(:,:,:),allocatable :: EUg,EUgder,CUg,& CUgder,DUg,Dugder,DtUg2,DtUg2der !(2,2,maxres) ! This common block contains vectors and matrices dependent on two @@ -87,6 +87,7 @@ real(kind=8),dimension(:,:,:,:),allocatable :: Ug2DtEUgder,& DtUg2EUgder !(2,2,2,maxres) ! common /rotat_old/ + real(kind=8),dimension(4) :: gmuij,gmuij1,gmuij2,gmuji1,gmuji2 real(kind=8),dimension(:),allocatable :: costab,sintab,& costab2,sintab2 !(maxres) ! This common block contains dipole-interaction matrices and their @@ -246,7 +247,7 @@ ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,& ecorr3_nucl ! energies for ions - real(kind=8) :: ecation_prot,ecationcation + real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber ! energies for protein nucleic acid interaction real(kind=8) :: escbase,epepbase,escpho,epeppho @@ -293,7 +294,7 @@ ! print *,"Processor",myrank," BROADCAST iorder" ! FG master sets up the WEIGHTS_ array which will be broadcast to the ! FG slaves as WEIGHTS array. - ! weights_(1)=wsc + weights_(1)=wsc weights_(2)=wscp weights_(3)=welec weights_(4)=wcorr @@ -328,8 +329,9 @@ weights_(41)=wcatcat weights_(42)=wcatprot weights_(46)=wscbase - weights_(47)=wscpho - weights_(48)=wpeppho + weights_(47)=wpepbase + weights_(48)=wscpho + weights_(49)=wpeppho ! wcatcat= weights(41) ! wcatprot=weights(42) @@ -375,8 +377,16 @@ wcatcat= weights(41) wcatprot=weights(42) wscbase=weights(46) - wscpho=weights(47) - wpeppho=weights(48) + wpepbase=weights(47) + wscpho=weights(48) + wpeppho=weights(49) +! welpsb=weights(28)*fact(1) +! +! wcorr_nucl= weights(37)*fact(1) +! wcorr3_nucl=weights(38)*fact(2) +! wtor_nucl= weights(35)*fact(1) +! wtor_d_nucl=weights(36)*fact(2) + endif time_Bcast=time_Bcast+MPI_Wtime()-time00 time_Bcastw=time_Bcastw+MPI_Wtime()-time00 @@ -609,9 +619,9 @@ .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0 & .or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then #endif - write(iout,*),"just befor eelec call" +! print *,"just befor eelec call" call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) -! write (iout,*) "ELEC calc" +! print *, "ELEC calc" else ees=0.0d0 evdw1=0.0d0 @@ -662,12 +672,26 @@ ! Calculate the virtual-bond-angle energy. ! write(iout,*) "in etotal afer edis",ipot - if (wang.gt.0.0d0) then - call ebend(ebe,ethetacnstr) +! if (wang.gt.0.0d0) then +! call ebend(ebe,ethetacnstr) +! else +! ebe=0 +! ethetacnstr=0 +! endif + if (wang.gt.0d0) then + if (tor_mode.eq.0) then + call ebend(ebe) + else +!C ebend kcc is Kubo cumulant clustered rigorous attemp to derive the +!C energy function + call ebend_kcc(ebe) + endif else - ebe=0 - ethetacnstr=0 + ebe=0.0d0 endif + ethetacnstr=0.0d0 + if (with_theta_constr) call etheta_constr(ethetacnstr) + ! write(iout,*) "in etotal afer ebe",ipot ! print *,"Processor",myrank," computed UB" @@ -681,12 +705,27 @@ ! Calculate the virtual-bond torsional energy. ! !d print *,'nterm=',nterm - if (wtor.gt.0) then - call etor(etors,edihcnstr) +! if (wtor.gt.0) then +! call etor(etors,edihcnstr) +! else +! etors=0 +! edihcnstr=0 +! endif + if (wtor.gt.0.0d0) then + if (tor_mode.eq.0) then + call etor(etors) + else +!C etor kcc is Kubo cumulant clustered rigorous attemp to derive the +!C energy function + call etor_kcc(etors) + endif else - etors=0 - edihcnstr=0 + etors=0.0d0 endif + edihcnstr=0.0d0 + if (ndih_constr.gt.0) call etor_constr(edihcnstr) +!c print *,"Processor",myrank," computed Utor" + ! print *,"Processor",myrank," computed Utor" ! @@ -763,6 +802,8 @@ call AFMforce(Eafmforce) else if (selfguide.gt.0) then call AFMvel(Eafmforce) + else + Eafmforce=0.0d0 endif endif if (tubemode.eq.1) then @@ -801,6 +842,7 @@ eespp=0.0d0 endif ! write(iout,*) ecorr_nucl,"ecorr_nucl",nres_molec(2) +! print *,"before ecatcat",wcatcat if (nfgtasks.gt.1) then if (fg_rank.eq.0) then call ecatcat(ecationcation) @@ -809,6 +851,7 @@ call ecatcat(ecationcation) endif call ecat_prot(ecation_prot) + call ecats_prot_amber(ecations_prot_amber) if (nres_molec(2).gt.0) then call eprot_sc_base(escbase) call epep_sc_base(epepbase) @@ -821,7 +864,7 @@ epeppho=0.0 endif ! call ecatcat(ecationcation) -! print *,"after ebend", ebe_nucl +! print *,"after ebend", wtor_nucl #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 #endif @@ -885,12 +928,13 @@ ! 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(42)=ecation_prot + energia(41)=ecationcation energia(46)=escbase energia(47)=epepbase energia(48)=escpho energia(49)=epeppho + energia(50)=ecations_prot_amber call sum_energy(energia,.true.) if (dyn_ss) call dyn_set_nss ! print *," Processor",myrank," left SUM_ENERGY" @@ -933,7 +977,7 @@ 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) :: ecation_prot,ecationcation,ecations_prot_amber real(kind=8) :: escbase,epepbase,escpho,epeppho integer :: i #ifdef MPI @@ -1011,12 +1055,14 @@ etors_d_nucl=energia(36) ecorr_nucl=energia(37) ecorr3_nucl=energia(38) - ecation_prot=energia(41) - ecationcation=energia(42) + ecation_prot=energia(42) + ecationcation=energia(41) escbase=energia(46) epepbase=energia(47) escpho=energia(48) epeppho=energia(49) + ecations_prot_amber=energia(50) + ! energia(41)=ecation_prot ! energia(42)=ecationcation @@ -1034,7 +1080,7 @@ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl& +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl& +wcatprot*ecation_prot+wcatcat*ecationcation+wscbase*escbase& - +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho + +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ecations_prot_amber #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc & @@ -1048,7 +1094,7 @@ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl& +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl& +wcatprot*ecation_prot+wcatcat*ecationcation+wscbase*escbase& - +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho + +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ecations_prot_amber #endif energia(0)=etot ! detecting NaNQ @@ -1156,7 +1202,12 @@ wtor=weights(13)*fact(1) wtor_d=weights(14)*fact(2) wsccor=weights(21)*fact(1) - + welpsb=weights(28)*fact(1) + wcorr_nucl= weights(37)*fact(1) + wcorr3_nucl=weights(38)*fact(2) + wtor_nucl= weights(35)*fact(1) + wtor_d_nucl=weights(36)*fact(2) + wpepbase=weights(47)*fact(1) return end subroutine rescale_weights !----------------------------------------------------------------------------- @@ -1176,7 +1227,7 @@ 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) :: ecation_prot,ecationcation,ecations_prot_amber real(kind=8) :: escbase,epepbase,escpho,epeppho etot=energia(0) @@ -1224,12 +1275,13 @@ etors_d_nucl=energia(36) ecorr_nucl=energia(37) ecorr3_nucl=energia(38) - ecation_prot=energia(41) - ecationcation=energia(42) + ecation_prot=energia(42) + ecationcation=energia(41) escbase=energia(46) epepbase=energia(47) escpho=energia(48) epeppho=energia(49) + ecations_prot_amber=energia(50) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,& estr,wbond,ebe,wang,& @@ -1245,7 +1297,7 @@ etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,& ecorr3_nucl,wcorr3_nucl,ecation_prot,wcatprot,ecationcation,wcatcat, & escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,& - etot + ecations_prot_amber,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ & @@ -1300,15 +1352,15 @@ ecorr,wcorr,& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,& - ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, & + ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforce, & etube,wtube, & estr_nucl,wbond_nucl, ebe_nucl,wang_nucl,& - evdwpp,wvdwpp_nucl,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb& - evdwsb,wvdwsb,eelsb,welsb,esbloc,wsbloc,etors_nucl,wtor_nucl& + 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,ecation_prot,wcatprot,ecationcation,wcatcat, & escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,& - etot + ecations_prot_amber,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ & @@ -1982,8 +2034,8 @@ ! write(iout,*)"c ", c(1,:), c(2,:), c(3,:) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) - sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj))) - sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj))) + sss_ele_cut=sscale_ele(1.0d0/(rij)) + sss_ele_grad=sscagrad_ele(1.0d0/(rij)) ! print *,sss_ele_cut,sss_ele_grad,& ! 1.0d0/(rij),r_cut_ele,rlamb_ele if (sss_ele_cut.le.0.0) cycle @@ -2045,7 +2097,7 @@ fac=rij*fac ! print *,'before fac',fac,rij,evdwij fac=fac+evdwij*sss_ele_grad/sss_ele_cut& - /sigma(itypi,itypj)*rij + *rij ! print *,'grad part scale',fac, & ! evdwij*sss_ele_grad/sss_ele_cut & ! /sigma(itypi,itypj)*rij @@ -2689,18 +2741,174 @@ ! include 'COMMON.FFIELD' real(kind=8) :: auxvec(2),auxmat(2,2) integer :: i,iti1,iti,k,l - real(kind=8) :: sin1,cos1,sin2,cos2,dwacos2,dwasin2 + real(kind=8) :: sin1,cos1,sin2,cos2,dwacos2,dwasin2,cost1,sint1,& + sint1sq,sint1cub,sint1cost1,b1k,b2k,aux ! print *,"in set matrices" ! ! Compute the virtual-bond-torsional-angle dependent quantities needed ! to calculate the el-loc multibody terms of various order. ! !AL el mu=0.0d0 + +#ifdef PARMAT + do i=ivec_start+2,ivec_end+2 +#else + do i=3,nres+1 +#endif + if (i.gt. nnt+2 .and. i.lt.nct+2) then + if (itype(i-2,1).eq.0) then + iti = nloctyp + else + iti = itype2loc(itype(i-2,1)) + endif + else + iti=nloctyp + endif +!c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then + if (i.gt. nnt+1 .and. i.lt.nct+1) then + iti1 = itype2loc(itype(i-1,1)) + else + iti1=nloctyp + endif +! print *,i,itype(i-2,1),iti +#ifdef NEWCORR + cost1=dcos(theta(i-1)) + sint1=dsin(theta(i-1)) + sint1sq=sint1*sint1 + sint1cub=sint1sq*sint1 + sint1cost1=2*sint1*cost1 +! print *,"cost1",cost1,theta(i-1) +!c write (iout,*) "bnew1",i,iti +!c write (iout,*) (bnew1(k,1,iti),k=1,3) +!c write (iout,*) (bnew1(k,2,iti),k=1,3) +!c write (iout,*) "bnew2",i,iti +!c write (iout,*) (bnew2(k,1,iti),k=1,3) +!c write (iout,*) (bnew2(k,2,iti),k=1,3) + k=1 +! print *,bnew1(1,k,iti),"bnew1" + do k=1,2 + b1k=bnew1(1,k,iti)+(bnew1(2,k,iti)+bnew1(3,k,iti)*cost1)*cost1 +! print *,b1k +! write(*,*) shape(b1) +! if(.not.allocated(b1)) print *, "WTF?" + b1(k,i-2)=sint1*b1k +! +! print *,b1(k,i-2) + + gtb1(k,i-2)=cost1*b1k-sint1sq*& + (bnew1(2,k,iti)+2*bnew1(3,k,iti)*cost1) +! print *,gtb1(k,i-2) + + b2k=bnew2(1,k,iti)+(bnew2(2,k,iti)+bnew2(3,k,iti)*cost1)*cost1 + b2(k,i-2)=sint1*b2k +! print *,b2(k,i-2) + + gtb2(k,i-2)=cost1*b2k-sint1sq*& + (bnew2(2,k,iti)+2*bnew2(3,k,iti)*cost1) +! print *,gtb2(k,i-2) + + enddo +! print *,b1k,b2k + do k=1,2 + aux=ccnew(1,k,iti)+(ccnew(2,k,iti)+ccnew(3,k,iti)*cost1)*cost1 + cc(1,k,i-2)=sint1sq*aux + gtcc(1,k,i-2)=sint1cost1*aux-sint1cub*& + (ccnew(2,k,iti)+2*ccnew(3,k,iti)*cost1) + aux=ddnew(1,k,iti)+(ddnew(2,k,iti)+ddnew(3,k,iti)*cost1)*cost1 + dd(1,k,i-2)=sint1sq*aux + gtdd(1,k,i-2)=sint1cost1*aux-sint1cub*& + (ddnew(2,k,iti)+2*ddnew(3,k,iti)*cost1) + enddo +! print *,"after cc" + cc(2,1,i-2)=cc(1,2,i-2) + cc(2,2,i-2)=-cc(1,1,i-2) + gtcc(2,1,i-2)=gtcc(1,2,i-2) + gtcc(2,2,i-2)=-gtcc(1,1,i-2) + dd(2,1,i-2)=dd(1,2,i-2) + dd(2,2,i-2)=-dd(1,1,i-2) + gtdd(2,1,i-2)=gtdd(1,2,i-2) + gtdd(2,2,i-2)=-gtdd(1,1,i-2) +! print *,"after dd" + + do k=1,2 + do l=1,2 + aux=eenew(1,l,k,iti)+eenew(2,l,k,iti)*cost1 + EE(l,k,i-2)=sint1sq*aux + gtEE(l,k,i-2)=sint1cost1*aux-sint1cub*eenew(2,l,k,iti) + enddo + enddo + EE(1,1,i-2)=EE(1,1,i-2)+e0new(1,iti)*cost1 + EE(1,2,i-2)=EE(1,2,i-2)+e0new(2,iti)+e0new(3,iti)*cost1 + EE(2,1,i-2)=EE(2,1,i-2)+e0new(2,iti)*cost1+e0new(3,iti) + EE(2,2,i-2)=EE(2,2,i-2)-e0new(1,iti) + gtEE(1,1,i-2)=gtEE(1,1,i-2)-e0new(1,iti)*sint1 + gtEE(1,2,i-2)=gtEE(1,2,i-2)-e0new(3,iti)*sint1 + gtEE(2,1,i-2)=gtEE(2,1,i-2)-e0new(2,iti)*sint1 +! print *,"after ee" + +!c b1tilde(1,i-2)=b1(1,i-2) +!c b1tilde(2,i-2)=-b1(2,i-2) +!c b2tilde(1,i-2)=b2(1,i-2) +!c b2tilde(2,i-2)=-b2(2,i-2) +#ifdef DEBUG + write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2) + write(iout,*) 'b1=',(b1(k,i-2),k=1,2) + write(iout,*) 'b2=',(b2(k,i-2),k=1,2) + write (iout,*) 'theta=', theta(i-1) +#endif +#else + if (i.gt. nnt+2 .and. i.lt.nct+2) then +! write(iout,*) "i,",molnum(i) +! print *, "i,",molnum(i),i,itype(i-2,1) + if (molnum(i).eq.1) then + iti = itype2loc(itype(i-2,1)) + else + iti=nloctyp + endif + else + iti=nloctyp + endif +!c write (iout,*) "i",i-1," itype",itype(i-2)," iti",iti +!c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then + if (i.gt. nnt+1 .and. i.lt.nct+1) then + iti1 = itype2loc(itype(i-1,1)) + else + iti1=nloctyp + endif +! print *,i,iti + b1(1,i-2)=b(3,iti) + b1(2,i-2)=b(5,iti) + b2(1,i-2)=b(2,iti) + b2(2,i-2)=b(4,iti) + do k=1,2 + do l=1,2 + CC(k,l,i-2)=ccold(k,l,iti) + DD(k,l,i-2)=ddold(k,l,iti) + EE(k,l,i-2)=eeold(k,l,iti) + enddo + enddo +#endif + b1tilde(1,i-2)= b1(1,i-2) + b1tilde(2,i-2)=-b1(2,i-2) + b2tilde(1,i-2)= b2(1,i-2) + b2tilde(2,i-2)=-b2(2,i-2) +!c + Ctilde(1,1,i-2)= CC(1,1,i-2) + Ctilde(1,2,i-2)= CC(1,2,i-2) + Ctilde(2,1,i-2)=-CC(2,1,i-2) + Ctilde(2,2,i-2)=-CC(2,2,i-2) +!c + Dtilde(1,1,i-2)= DD(1,1,i-2) + Dtilde(1,2,i-2)= DD(1,2,i-2) + Dtilde(2,1,i-2)=-DD(2,1,i-2) + Dtilde(2,2,i-2)=-DD(2,2,i-2) + enddo #ifdef PARMAT do i=ivec_start+2,ivec_end+2 #else do i=3,nres+1 #endif + ! print *,i,"i" if (i .lt. nres+1) then sin1=dsin(phi(i)) @@ -2773,37 +2981,43 @@ if (itype(i-2,1).eq.0) then iti=ntortyp+1 else - iti = itortyp(itype(i-2,1)) + iti = itype2loc(itype(i-2,1)) endif else - iti=ntortyp+1 + iti=nloctyp 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 + iti1=nloctyp else - iti1 = itortyp(itype(i-1,1)) + iti1 = itype2loc(itype(i-1,1)) endif else - iti1=ntortyp+1 + iti1=nloctyp endif ! print *,iti,i,"iti",iti1,itype(i-1,1),itype(i-2,1) !d write (iout,*) '*******i',i,' iti1',iti -!d write (iout,*) 'b1',b1(:,iti) -!d write (iout,*) 'b2',b2(:,iti) +! write (iout,*) 'b1',b1(:,iti) +! write (iout,*) 'b2',b2(:,i-2) !d write (iout,*) 'Ug',Ug(:,:,i-2) ! if (i .gt. iatel_s+2) then if (i .gt. nnt+2) then - call matvec2(Ug(1,1,i-2),b2(1,iti),Ub2(1,i-2)) - call matmat2(EE(1,1,iti),Ug(1,1,i-2),EUg(1,1,i-2)) + call matvec2(Ug(1,1,i-2),b2(1,i-2),Ub2(1,i-2)) +#ifdef NEWCORR + call matvec2(Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2)) +!c write (iout,*) Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2),"chuj" +#endif + + call matmat2(EE(1,1,i-2),Ug(1,1,i-2),EUg(1,1,i-2)) + call matmat2(gtEE(1,1,i-2),Ug(1,1,i-2),gtEUg(1,1,i-2)) if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) & then - call matmat2(CC(1,1,iti),Ug(1,1,i-2),CUg(1,1,i-2)) - call matmat2(DD(1,1,iti),Ug(1,1,i-2),DUg(1,1,i-2)) - call matmat2(Dtilde(1,1,iti),Ug2(1,1,i-2),DtUg2(1,1,i-2)) - call matvec2(Ctilde(1,1,iti1),obrot(1,i-2),Ctobr(1,i-2)) - call matvec2(Dtilde(1,1,iti),obrot2(1,i-2),Dtobr2(1,i-2)) + call matmat2(CC(1,1,i-2),Ug(1,1,i-2),CUg(1,1,i-2)) + call matmat2(DD(1,1,i-2),Ug(1,1,i-2),DUg(1,1,i-2)) + call matmat2(Dtilde(1,1,i-2),Ug2(1,1,i-2),DtUg2(1,1,i-2)) + call matvec2(Ctilde(1,1,i-1),obrot(1,i-2),Ctobr(1,i-2)) + call matvec2(Dtilde(1,1,i-2),obrot2(1,i-2),Dtobr2(1,i-2)) endif else do k=1,2 @@ -2818,44 +3032,44 @@ enddo enddo endif - call matvec2(Ugder(1,1,i-2),b2(1,iti),Ub2der(1,i-2)) - call matmat2(EE(1,1,iti),Ugder(1,1,i-2),EUgder(1,1,i-2)) + call matvec2(Ugder(1,1,i-2),b2(1,i-2),Ub2der(1,i-2)) + call matmat2(EE(1,1,i-2),Ugder(1,1,i-2),EUgder(1,1,i-2)) do k=1,2 muder(k,i-2)=Ub2der(k,i-2) 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).eq.0) then - iti1=ntortyp+1 + iti1=nloctyp elseif (itype(i-1,1).le.ntyp) then - iti1 = itortyp(itype(i-1,1)) + iti1 = itype2loc(itype(i-1,1)) else - iti1=ntortyp+1 + iti1=nloctyp endif else - iti1=ntortyp+1 + iti1=nloctyp endif do k=1,2 - mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1) + mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1) enddo -! if (energy_dec) write (iout,*) 'Ub2 ',i,Ub2(:,i-2) -! if (energy_dec) write (iout,*) 'b1 ',iti1,b1(:,iti1) -! if (energy_dec) write (iout,*) 'mu ',i,iti1,mu(:,i-2) + if (energy_dec) write (iout,*) 'Ub2 ',i,Ub2(:,i-2) + if (energy_dec) write (iout,*) 'b1 ',iti1,b1(:,i-1) + if (energy_dec) write (iout,*) 'mu ',i,iti1,mu(:,i-2) !d write (iout,*) 'mu1',mu1(:,i-2) !d write (iout,*) 'mu2',mu2(:,i-2) if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) & then - call matmat2(CC(1,1,iti1),Ugder(1,1,i-2),CUgder(1,1,i-2)) - call matmat2(DD(1,1,iti),Ugder(1,1,i-2),DUgder(1,1,i-2)) - call matmat2(Dtilde(1,1,iti),Ug2der(1,1,i-2),DtUg2der(1,1,i-2)) - call matvec2(Ctilde(1,1,iti1),obrot_der(1,i-2),Ctobrder(1,i-2)) - call matvec2(Dtilde(1,1,iti),obrot2_der(1,i-2),Dtobr2der(1,i-2)) + call matmat2(CC(1,1,i-1),Ugder(1,1,i-2),CUgder(1,1,i-2)) + call matmat2(DD(1,1,i-2),Ugder(1,1,i-2),DUgder(1,1,i-2)) + call matmat2(Dtilde(1,1,i-2),Ug2der(1,1,i-2),DtUg2der(1,1,i-2)) + call matvec2(Ctilde(1,1,i-1),obrot_der(1,i-2),Ctobrder(1,i-2)) + call matvec2(Dtilde(1,1,i-2),obrot2_der(1,i-2),Dtobr2der(1,i-2)) ! Vectors and matrices dependent on a single virtual-bond dihedral. - call matvec2(DD(1,1,iti),b1tilde(1,iti1),auxvec(1)) + call matvec2(DD(1,1,i-2),b1tilde(1,i-1),auxvec(1)) call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2)) call matvec2(Ug2der(1,1,i-2),auxvec(1),Ug2Db1tder(1,i-2)) - call matvec2(CC(1,1,iti1),Ub2(1,i-2),CUgb2(1,i-2)) - call matvec2(CC(1,1,iti1),Ub2der(1,i-2),CUgb2der(1,i-2)) + call matvec2(CC(1,1,i-1),Ub2(1,i-2),CUgb2(1,i-2)) + call matvec2(CC(1,1,i-1),Ub2der(1,i-2),CUgb2der(1,i-2)) call matmat2(EUg(1,1,i-2),CC(1,1,iti1),EUgC(1,1,i-2)) call matmat2(EUgder(1,1,i-2),CC(1,1,iti1),EUgCder(1,1,i-2)) call matmat2(EUg(1,1,i-2),DD(1,1,iti1),EUgD(1,1,i-2)) @@ -3476,6 +3690,7 @@ real(kind=8),dimension(2,2) :: acipa !el,a_temp !el real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1 real(kind=8),dimension(4) :: muij + real(kind=8) :: geel_loc_ij,geel_loc_ji real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,& dist_temp, dist_init,rlocshield,fracinbuf integer xshift,yshift,zshift,ilist,iresshield @@ -3917,6 +4132,15 @@ do l=1,2 kkk=kkk+1 muij(kkk)=mu(k,i)*mu(l,j) +#ifdef NEWCORR + gmuij1(kkk)=gtb1(k,i+1)*mu(l,j) +!c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j) + gmuij2(kkk)=gUb2(k,i)*mu(l,j) + gmuji1(kkk)=mu(k,i)*gtb1(l,j+1) +!c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i) + gmuji2(kkk)=mu(k,i)*gUb2(l,j) +#endif + enddo enddo !d write (iout,*) 'EELEC: i',i,' j',j @@ -4143,6 +4367,61 @@ enddo endif +#ifdef NEWCORR + geel_loc_ij=(a22*gmuij1(1)& + +a23*gmuij1(2)& + +a32*gmuij1(3)& + +a33*gmuij1(4))& + *fac_shield(i)*fac_shield(j)& + *sss_ele_cut + +!c write(iout,*) "derivative over thatai" +!c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), +!c & a33*gmuij1(4) + gloc(nphi+i,icg)=gloc(nphi+i,icg)+& + geel_loc_ij*wel_loc +!c write(iout,*) "derivative over thatai-1" +!c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3), +!c & a33*gmuij2(4) + geel_loc_ij=& + a22*gmuij2(1)& + +a23*gmuij2(2)& + +a32*gmuij2(3)& + +a33*gmuij2(4) + gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+& + geel_loc_ij*wel_loc& + *fac_shield(i)*fac_shield(j)& + *sss_ele_cut + + +!c Derivative over j residue + geel_loc_ji=a22*gmuji1(1)& + +a23*gmuji1(2)& + +a32*gmuji1(3)& + +a33*gmuji1(4) +!c write(iout,*) "derivative over thataj" +!c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3), +!c & a33*gmuji1(4) + + gloc(nphi+j,icg)=gloc(nphi+j,icg)+& + geel_loc_ji*wel_loc& + *fac_shield(i)*fac_shield(j)& + *sss_ele_cut + + + geel_loc_ji=& + +a22*gmuji2(1)& + +a23*gmuji2(2)& + +a32*gmuji2(3)& + +a33*gmuji2(4) +!c write(iout,*) "derivative over thataj-1" +!c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), +!c & a33*gmuji2(4) + gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+& + geel_loc_ji*wel_loc& + *fac_shield(i)*fac_shield(j)& + *sss_ele_cut +#endif ! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij ! eel_loc_ij=0.0 @@ -4501,7 +4780,9 @@ ! include 'COMMON.CONTROL' real(kind=8),dimension(3) :: ggg real(kind=8),dimension(2,2) :: auxmat,auxmat1,auxmat2,pizda,& - e1t,e2t,e3t,e1tder,e2tder,e3tder,e1a,ae3,ae3e2 + e1t,e2t,e3t,e1tder,e2tder,e3tder,e1a,ae3,ae3e2,gpizda1,& + gpizda2,auxgmat1,auxgmatt1,auxgmat2,auxgmatt2 + real(kind=8),dimension(2) :: auxvec,auxvec1 !el real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1 real(kind=8),dimension(2,2) :: auxmat3 !el, a_temp @@ -4561,8 +4842,15 @@ !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC !d call checkint_turn3(i,a_temp,eello_turn3_num) call matmat2(EUg(1,1,i+1),EUg(1,1,i+2),auxmat(1,1)) + call matmat2(gtEUg(1,1,i+1),EUg(1,1,i+2),auxgmat1(1,1)) + call matmat2(EUg(1,1,i+1),gtEUg(1,1,i+2),auxgmat2(1,1)) call transpose2(auxmat(1,1),auxmat1(1,1)) + call transpose2(auxgmat1(1,1),auxgmatt1(1,1)) + call transpose2(auxgmat2(1,1),auxgmatt2(1,1)) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) + call matmat2(a_temp(1,1),auxgmatt1(1,1),gpizda1(1,1)) + call matmat2(a_temp(1,1),auxgmatt2(1,1),gpizda2(1,1)) + if (shield_mode.eq.0) then fac_shield(i)=1.0d0 fac_shield(j)=1.0d0 @@ -4577,6 +4865,18 @@ if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2)) +!C#ifdef NEWCORR +!C Derivatives in theta + gloc(nphi+i,icg)=gloc(nphi+i,icg) & + +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3& + *fac_shield(i)*fac_shield(j) + gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)& + +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3& + *fac_shield(i)*fac_shield(j) +!C#endif + + + if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & (shield_mode.gt.0)) then !C print *,i,j @@ -4714,8 +5014,14 @@ ! include 'COMMON.CONTROL' real(kind=8),dimension(3) :: ggg real(kind=8),dimension(2,2) :: auxmat,auxmat1,auxmat2,pizda,& - e1t,e2t,e3t,e1tder,e2tder,e3tder,e1a,ae3,ae3e2 - real(kind=8),dimension(2) :: auxvec,auxvec1 + e1t,e2t,e3t,e1tder,e2tder,e3tder,e1a,ae3,ae3e2,& + gte1t,gte2t,gte3t,& + gte1a,gtae3,gtae3e2, ae3gte2,& + gtEpizda1,gtEpizda2,gtEpizda3 + + real(kind=8),dimension(2) :: auxvec,auxvec1,auxgEvec1,auxgEvec2,& + auxgEvec3,auxgvec + !el real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1 real(kind=8),dimension(2,2) :: auxmat3 !el a_temp !el real(kind=8) :: a22,a23,a32,a33,dxi,dyi,dzi,dx_normi,dy_normi,& @@ -4727,7 +5033,7 @@ !el local variables integer :: i,j,iti1,iti2,iti3,l,k,ilist,iresshield real(kind=8) :: eello_turn4,s1,s2,s3,zj,fracinbuf,eello_t4,& - rlocshield + rlocshield,gs23,gs32,gsE13,gs13,gs21,gsE31,gsEE1,gsEE2,gsEE3 j=i+3 ! if (j.ne.20) return @@ -4775,22 +5081,63 @@ a_temp(1,2)=a23 a_temp(2,1)=a32 a_temp(2,2)=a33 - iti1=itortyp(itype(i+1,1)) - iti2=itortyp(itype(i+2,1)) - iti3=itortyp(itype(i+3,1)) + iti1=i+1 + iti2=i+2 + iti3=i+3 ! write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3 call transpose2(EUg(1,1,i+1),e1t(1,1)) call transpose2(Eug(1,1,i+2),e2t(1,1)) call transpose2(Eug(1,1,i+3),e3t(1,1)) +!C Ematrix derivative in theta + call transpose2(gtEUg(1,1,i+1),gte1t(1,1)) + call transpose2(gtEug(1,1,i+2),gte2t(1,1)) + call transpose2(gtEug(1,1,i+3),gte3t(1,1)) + call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1)) call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1)) + call matmat2(gte1t(1,1),a_temp(1,1),gte1a(1,1)) + call matvec2(e1a(1,1),gUb2(1,i+3),auxgvec(1)) +!c auxalary matrix of E i+1 + call matvec2(gte1a(1,1),Ub2(1,i+3),auxgEvec1(1)) s1=scalar2(b1(1,iti2),auxvec(1)) +!c derivative of theta i+2 with constant i+3 + gs23=scalar2(gtb1(1,i+2),auxvec(1)) +!c derivative of theta i+2 with constant i+2 + gs32=scalar2(b1(1,i+2),auxgvec(1)) +!c derivative of E matix in theta of i+1 + gsE13=scalar2(b1(1,i+2),auxgEvec1(1)) + call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1)) + call matmat2(a_temp(1,1),gte3t(1,1),gtae3(1,1)) call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) - s2=scalar2(b1(1,iti1),auxvec(1)) +!c auxilary matrix auxgvec of Ub2 with constant E matirx + call matvec2(ae3(1,1),gUb2(1,i+2),auxgvec(1)) +!c auxilary matrix auxgEvec1 of E matix with Ub2 constant + call matvec2(gtae3(1,1),Ub2(1,i+2),auxgEvec3(1)) + s2=scalar2(b1(1,i+1),auxvec(1)) +!c derivative of theta i+1 with constant i+3 + gs13=scalar2(gtb1(1,i+1),auxvec(1)) +!c derivative of theta i+2 with constant i+1 + gs21=scalar2(b1(1,i+1),auxgvec(1)) +!c derivative of theta i+3 with constant i+1 + gsE31=scalar2(b1(1,i+1),auxgEvec3(1)) + call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) + call matmat2(gtae3(1,1),e2t(1,1),gtae3e2(1,1)) +!c ae3gte2 is derivative over i+2 + call matmat2(ae3(1,1),gte2t(1,1),ae3gte2(1,1)) + call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) + call matmat2(ae3e2(1,1),gte1t(1,1),gtEpizda1(1,1)) +!c i+2 + call matmat2(ae3gte2(1,1),e1t(1,1),gtEpizda2(1,1)) +!c i+3 + call matmat2(gtae3e2(1,1),e1t(1,1),gtEpizda3(1,1)) + s3=0.5d0*(pizda(1,1)+pizda(2,2)) + gsEE1=0.5d0*(gtEpizda1(1,1)+gtEpizda1(2,2)) + gsEE2=0.5d0*(gtEpizda2(1,1)+gtEpizda2(2,2)) + gsEE3=0.5d0*(gtEpizda3(1,1)+gtEpizda3(2,2)) if (shield_mode.eq.0) then fac_shield(i)=1.0 fac_shield(j)=1.0 @@ -4844,7 +5191,21 @@ ! print *,"gshieldc_t4(k,j+1)",j,gshieldc_t4(k,j+1) enddo endif - +#ifdef NEWCORR + gloc(nphi+i,icg)=gloc(nphi+i,icg)& + -(gs13+gsE13+gsEE1)*wturn4& + *fac_shield(i)*fac_shield(j) + gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg)& + -(gs23+gs21+gsEE2)*wturn4& + *fac_shield(i)*fac_shield(j) + + gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg)& + -(gs32+gsE31+gsEE3)*wturn4& + *fac_shield(i)*fac_shield(j) + +!c gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)- +!c & gs2 +#endif if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eturn4',i,j,-(s1+s2+s3) !d write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3), @@ -4853,7 +5214,7 @@ call transpose2(EUgder(1,1,i+1),e1tder(1,1)) call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1)) call matvec2(auxmat(1,1),Ub2(1,i+3),auxvec(1)) - s1=scalar2(b1(1,iti2),auxvec(1)) + s1=scalar2(b1(1,i+1),auxvec(1)) call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3) & @@ -5897,7 +6258,7 @@ end subroutine theteng #else !----------------------------------------------------------------------------- - subroutine ebend(etheta,ethetacnstr) + subroutine ebend(etheta) ! ! Evaluate the virtual-bond-angle energy given the virtual-bond dihedral ! angles gamma and its derivatives in consecutive thetas and gammas. @@ -6103,34 +6464,6 @@ enddo !-----------thete constrains ! if (tor_mode.ne.2) then - ethetacnstr=0.0d0 - print *,ithetaconstr_start,ithetaconstr_end,"TU" - do i=ithetaconstr_start,ithetaconstr_end - itheta=itheta_constr(i) - thetiii=theta(itheta) - difi=pinorm(thetiii-theta_constr0(i)) - if (difi.gt.theta_drange(i)) then - difi=difi-theta_drange(i) - ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 - gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) & - +for_thet_constr(i)*difi**3 - else if (difi.lt.-drange(i)) then - difi=difi+drange(i) - ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 - gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) & - +for_thet_constr(i)*difi**3 - else - difi=0.0 - endif - if (energy_dec) then - write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", & - i,itheta,rad2deg*thetiii, & - rad2deg*theta_constr0(i), rad2deg*theta_drange(i), & - rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, & - gloc(itheta+nphi-2,icg) - endif - enddo -! endif return end subroutine ebend @@ -7019,7 +7352,7 @@ end subroutine etor_d #else !----------------------------------------------------------------------------- - subroutine etor(etors,edihcnstr) + subroutine etor(etors) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.VAR' @@ -7101,8 +7434,156 @@ ! write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) enddo ! 6/20/98 - dihedral angle constraints - edihcnstr=0.0d0 -! do i=1,ndih_constr + return + end subroutine etor +!C The rigorous attempt to derive energy function +!------------------------------------------------------------------------------------------- + subroutine etor_kcc(etors) + double precision c1(0:maxval_kcc),c2(0:maxval_kcc) + real(kind=8) :: etors,glocig,glocit1,glocit2,sinthet1,& + sinthet2,costhet1,costhet2,sint1t2,sint1t2n,phii,sinphi,cosphi,& + sint1t2n1,sumvalc,gradvalct1,gradvalct2,sumvals,gradvalst1,& + gradvalst2,etori + logical lprn + integer :: i,j,itori,itori1,nval,k,l + + if (lprn) write (iout,*) "etor_kcc tor_mode",tor_mode + etors=0.0D0 + do i=iphi_start,iphi_end +!C ANY TWO ARE DUMMY ATOMS in row CYCLE +!c if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or. +!c & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)) .or. +!c & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle + if (itype(i-2,1).eq.ntyp1.or. itype(i-1,1).eq.ntyp1 & + .or. itype(i,1).eq.ntyp1 .or. itype(i-3,1).eq.ntyp1) cycle + itori=itortyp(itype(i-2,1)) + itori1=itortyp(itype(i-1,1)) + phii=phi(i) + glocig=0.0D0 + glocit1=0.0d0 + glocit2=0.0d0 +!C to avoid multiple devision by 2 +!c theti22=0.5d0*theta(i) +!C theta 12 is the theta_1 /2 +!C theta 22 is theta_2 /2 +!c theti12=0.5d0*theta(i-1) +!C and appropriate sinus function + sinthet1=dsin(theta(i-1)) + sinthet2=dsin(theta(i)) + costhet1=dcos(theta(i-1)) + costhet2=dcos(theta(i)) +!C to speed up lets store its mutliplication + sint1t2=sinthet2*sinthet1 + sint1t2n=1.0d0 +!C \sum_{i=1}^n (sin(theta_1) * sin(theta_2))^n * (c_n* cos(n*gamma) +!C +d_n*sin(n*gamma)) * +!C \sum_{i=1}^m (1+a_m*Tb_m(cos(theta_1 /2))+b_m*Tb_m(cos(theta_2 /2))) +!C we have two sum 1) Non-Chebyshev which is with n and gamma + nval=nterm_kcc_Tb(itori,itori1) + c1(0)=0.0d0 + c2(0)=0.0d0 + c1(1)=1.0d0 + c2(1)=1.0d0 + do j=2,nval + c1(j)=c1(j-1)*costhet1 + c2(j)=c2(j-1)*costhet2 + enddo + etori=0.0d0 + + do j=1,nterm_kcc(itori,itori1) + cosphi=dcos(j*phii) + sinphi=dsin(j*phii) + sint1t2n1=sint1t2n + sint1t2n=sint1t2n*sint1t2 + sumvalc=0.0d0 + gradvalct1=0.0d0 + gradvalct2=0.0d0 + do k=1,nval + do l=1,nval + sumvalc=sumvalc+v1_kcc(l,k,j,itori1,itori)*c1(k)*c2(l) + gradvalct1=gradvalct1+ & + (k-1)*v1_kcc(l,k,j,itori1,itori)*c1(k-1)*c2(l) + gradvalct2=gradvalct2+ & + (l-1)*v1_kcc(l,k,j,itori1,itori)*c1(k)*c2(l-1) + enddo + enddo + gradvalct1=-gradvalct1*sinthet1 + gradvalct2=-gradvalct2*sinthet2 + sumvals=0.0d0 + gradvalst1=0.0d0 + gradvalst2=0.0d0 + do k=1,nval + do l=1,nval + sumvals=sumvals+v2_kcc(l,k,j,itori1,itori)*c1(k)*c2(l) + gradvalst1=gradvalst1+ & + (k-1)*v2_kcc(l,k,j,itori1,itori)*c1(k-1)*c2(l) + gradvalst2=gradvalst2+ & + (l-1)*v2_kcc(l,k,j,itori1,itori)*c1(k)*c2(l-1) + enddo + enddo + gradvalst1=-gradvalst1*sinthet1 + gradvalst2=-gradvalst2*sinthet2 + if (lprn) write (iout,*)j,"sumvalc",sumvalc," sumvals",sumvals + etori=etori+sint1t2n*(sumvalc*cosphi+sumvals*sinphi) +!C glocig is the gradient local i site in gamma + glocig=glocig+j*sint1t2n*(sumvals*cosphi-sumvalc*sinphi) +!C now gradient over theta_1 + glocit1=glocit1+sint1t2n*(gradvalct1*cosphi+gradvalst1*sinphi)& + +j*sint1t2n1*costhet1*sinthet2*(sumvalc*cosphi+sumvals*sinphi) + glocit2=glocit2+sint1t2n*(gradvalct2*cosphi+gradvalst2*sinphi)& + +j*sint1t2n1*sinthet1*costhet2*(sumvalc*cosphi+sumvals*sinphi) + enddo ! j + etors=etors+etori + gloc(i-3,icg)=gloc(i-3,icg)+wtor*glocig +!C derivative over theta1 + gloc(nphi+i-3,icg)=gloc(nphi+i-3,icg)+wtor*glocit1 +!C now derivative over theta2 + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wtor*glocit2 + if (lprn) then + write (iout,*) i-2,i-1,itype(i-2,1),itype(i-1,1),itori,itori1,& + theta(i-1)*rad2deg,theta(i)*rad2deg,phii*rad2deg,etori + write (iout,*) "c1",(c1(k),k=0,nval), & + " c2",(c2(k),k=0,nval) + endif + enddo + return + end subroutine etor_kcc +!------------------------------------------------------------------------------ + + subroutine etor_constr(edihcnstr) + real(kind=8) :: etors,edihcnstr + logical :: lprn +!el local variables + integer :: i,j,iblock,itori,itori1 + real(kind=8) :: phii,gloci,v1ij,v2ij,cosphi,sinphi,& + vl1ij,vl2ij,vl3ij,pom1,difi,etors_ii,pom,& + gaudih_i,gauder_i,s,cos_i,dexpcos_i + + if (raw_psipred) then + do i=idihconstr_start,idihconstr_end + itori=idih_constr(i) + phii=phi(itori) + gaudih_i=vpsipred(1,i) + gauder_i=0.0d0 + do j=1,2 + s = sdihed(j,i) + cos_i=(1.0d0-dcos(phii-phibound(j,i)))/s**2 + dexpcos_i=dexp(-cos_i*cos_i) + gaudih_i=gaudih_i+vpsipred(j+1,i)*dexpcos_i + gauder_i=gauder_i-2*vpsipred(j+1,i)*dsin(phii-phibound(j,i)) & + *cos_i*dexpcos_i/s**2 + enddo + edihcnstr=edihcnstr-wdihc*dlog(gaudih_i) + gloc(itori-3,icg)=gloc(itori-3,icg)-wdihc*gauder_i/gaudih_i + if (energy_dec) & + write (iout,'(2i5,f8.3,f8.5,2(f8.5,2f8.3),f10.5)') & + i,itori,phii*rad2deg,vpsipred(1,i),vpsipred(2,i),& + phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,vpsipred(3,i),& + phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,& + -wdihc*dlog(gaudih_i) + enddo + else + do i=idihconstr_start,idihconstr_end itori=idih_constr(i) phii=phi(itori) @@ -7118,13 +7599,13 @@ else difi=0.0 endif -!d write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii, -!d & rad2deg*phi0(i), rad2deg*drange(i), -!d & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) enddo -!d write (iout,*) 'edihcnstr',edihcnstr + + endif + return - end subroutine etor + + end subroutine etor_constr !----------------------------------------------------------------------------- subroutine etor_d(etors_d) ! 6/23/01 Compute double torsional energy @@ -7215,6 +7696,75 @@ return end subroutine etor_d #endif + + subroutine ebend_kcc(etheta) + logical lprn + double precision thybt1(maxang_kcc),etheta + integer :: i,iti,j,ihelp + real (kind=8) :: sinthet,costhet,sumth1thyb,gradthybt1 +!C Set lprn=.true. for debugging + lprn=energy_dec +!c lprn=.true. +!C print *,"wchodze kcc" + if (lprn) write (iout,*) "ebend_kcc tor_mode",tor_mode + etheta=0.0D0 + do i=ithet_start,ithet_end +!c print *,i,itype(i-1),itype(i),itype(i-2) + if ((itype(i-1,1).eq.ntyp1).or.itype(i-2,1).eq.ntyp1 & + .or.itype(i,1).eq.ntyp1) cycle + iti=iabs(itortyp(itype(i-1,1))) + sinthet=dsin(theta(i)) + costhet=dcos(theta(i)) + do j=1,nbend_kcc_Tb(iti) + thybt1(j)=v1bend_chyb(j,iti) + enddo + sumth1thyb=v1bend_chyb(0,iti)+ & + tschebyshev(1,nbend_kcc_Tb(iti),thybt1(1),costhet) + if (lprn) write (iout,*) i-1,itype(i-1,1),iti,theta(i)*rad2deg,& + sumth1thyb + ihelp=nbend_kcc_Tb(iti)-1 + gradthybt1=gradtschebyshev(0,ihelp,thybt1(1),costhet) + etheta=etheta+sumth1thyb +!C print *,sumth1thyb,gradthybt1,sinthet*(-0.5d0) + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)-wang*gradthybt1*sinthet + enddo + return + end subroutine ebend_kcc +!c------------ +!c------------------------------------------------------------------------------------- + subroutine etheta_constr(ethetacnstr) + real (kind=8) :: ethetacnstr,thetiii,difi + integer :: i,itheta + ethetacnstr=0.0d0 +!C print *,ithetaconstr_start,ithetaconstr_end,"TU" + do i=ithetaconstr_start,ithetaconstr_end + itheta=itheta_constr(i) + thetiii=theta(itheta) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) & + +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) & + +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif + if (energy_dec) then + write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",& + i,itheta,rad2deg*thetiii,& + rad2deg*theta_constr0(i), rad2deg*theta_drange(i),& + rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,& + gloc(itheta+nphi-2,icg) + endif + enddo + return + end subroutine etheta_constr + !----------------------------------------------------------------------------- subroutine eback_sc_corr(esccor) ! 7/21/2007 Correlations between the backbone-local and side-chain-local @@ -7434,7 +7984,7 @@ jj,jp,kk,j1,jp1,jjc,iii,nnn,iproc ! Set lprn=.true. for debugging - lprn=.true. + lprn=.false. #ifdef MPI ! maxconts=nres/4 if(.not.allocated(zapas)) allocate(zapas(max_dim,maxconts,nfgtasks)) @@ -8411,9 +8961,9 @@ iti1 = itortyp(itype(i+1,1)) if (j.lt.nres-1) then - itj1 = itortyp(itype(j+1,1)) + itj1 = itype2loc(itype(j+1,1)) else - itj1=ntortyp+1 + itj1=nloctyp endif do iii=1,2 dipi(iii,1)=Ub2(iii,i) @@ -10685,6 +11235,7 @@ #ifdef TIMING time01=MPI_Wtime() #endif +!#define DEBUG #ifdef DEBUG write (iout,*) "sum_gradient gvdwc, gvdwx" do i=1,nres @@ -10718,18 +11269,18 @@ ! write (iout,'(i5,3f10.5,2x,f10.5)') ! & i,(gcorr4_turn(j,i),j=1,3),gel_loc_turn4(i) ! enddo - write (iout,*) "gvdwc gvdwc_scp gvdwc_scpp" - do i=1,nres - write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') & - i,(gvdwc(j,i),j=1,3),(gvdwc_scp(j,i),j=1,3),& - (gvdwc_scpp(j,i),j=1,3) - enddo - write (iout,*) "gelc_long gvdwpp gel_loc_long" - do i=1,nres - write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') & - i,(gelc_long(j,i),j=1,3),(gvdwpp(j,i),j=1,3),& - (gelc_loc_long(j,i),j=1,3) - enddo +! write (iout,*) "gvdwc gvdwc_scp gvdwc_scpp" +! do i=1,nres +! write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') & +! i,(gvdwc(j,i),j=1,3),(gvdwc_scp(j,i),j=1,3),& +! (gvdwc_scpp(j,i),j=1,3) +! enddo +! write (iout,*) "gelc_long gvdwpp gel_loc_long" +! do i=1,nres +! write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') & +! i,(gelc_long(j,i),j=1,3),(gvdwpp(j,i),j=1,3),& +! (gelc_loc_long(j,i),j=1,3) +! enddo call flush(iout) #endif #ifdef SPLITELE @@ -10797,7 +11348,7 @@ +wcorr3_nucl*gradcorr3_nucl(j,i) +& wcatprot* gradpepcat(j,i)+ & wcatcat*gradcatcat(j,i)+ & - wscbase*gvdwc_scbase(j,i) & + wscbase*gvdwc_scbase(j,i)+ & wpepbase*gvdwc_pepbase(j,i)+& wscpho*gvdwc_scpho(j,i)+& wpeppho*gvdwc_peppho(j,i) @@ -11028,7 +11579,7 @@ +gradafm(j,i) & +wliptran*gliptranc(j,i) & +welec*gshieldc(j,i) & - +welec*gshieldc_loc(j,) & + +welec*gshieldc_loc(j,i) & +wcorr*gshieldc_ec(j,i) & +wcorr*gshieldc_loc_ec(j,i) & +wturn3*gshieldc_t3(j,i) & @@ -11340,50 +11891,119 @@ enddo return end subroutine sc_grad -#ifdef CRYST_THETA -!----------------------------------------------------------------------------- - subroutine mixder(thetai,thet_pred_mean,theta0i,E_tc_t) - use comm_calcthet + subroutine sc_grad_cat ! implicit real*8 (a-h,o-z) + use calc_data ! include 'DIMENSIONS' -! include 'COMMON.LOCAL' +! include 'COMMON.CHAIN' +! include 'COMMON.DERIV' +! include 'COMMON.CALC' ! include 'COMMON.IOUNITS' -!el real(kind=8) :: term1,term2,termm,diffak,ratak,& -!el ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq,& -!el delthe0,sig0inv,sigtc,sigsqtc,delthec, - real(kind=8) :: thetai,thet_pred_mean,theta0i,E_tc_t - real(kind=8) :: t3,t6,t9,t12,t14,t16,t21,t23,t26,t27,t32,t40 -!el integer :: it -!el common /calcthet/ term1,term2,termm,diffak,ratak,& -!el ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq,& -!el delthe0,sig0inv,sigtc,sigsqtc,delthec,it -!el local variables + real(kind=8), dimension(3) :: dcosom1,dcosom2 +! print *,"wchodze" + eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 & + +dCAVdOM1+ dGCLdOM1+ dPOLdOM1 + eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 & + +dCAVdOM2+ dGCLdOM2+ dPOLdOM2 - delthec=thetai-thet_pred_mean - delthe0=thetai-theta0i -! "Thank you" to MAPLE (probably spared one day of hand-differentiation). - t3 = thetai-thet_pred_mean - t6 = t3**2 - t9 = term1 - t12 = t3*sigcsq - t14 = t12+t6*sigsqtc - t16 = 1.0d0 - t21 = thetai-theta0i - t23 = t21**2 - t26 = term2 - t27 = t21*t26 - t32 = termexp - t40 = t32**2 - E_tc_t = -((sigcsq+2.D0*t3*sigsqtc)*t9-t14*sigcsq*t3*t16*t9 & - -aktc*sig0inv*t27)/t32+(t14*t9+aktc*t26)/t40 & - *(-t12*t9-ak*sig0inv*t27) - return - end subroutine mixder -#endif -!----------------------------------------------------------------------------- -! cartder.F -!----------------------------------------------------------------------------- + eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 & + -2.0D0*alf12*eps3der+sigder*sigsq_om12& + +dCAVdOM12+ dGCLdOM12 +! diagnostics only +! eom1=0.0d0 +! eom2=0.0d0 +! eom12=evdwij*eps1_om12 +! end diagnostics +! write (iout,*) "eps2der",eps2der," eps3der",eps3der,& +! " sigder",sigder +! write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12 +! write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12 +!C print *,sss_ele_cut,'in sc_grad' + + do k=1,3 + dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k)) + dcosom2(k)=rij*(dc_norm(k,j)-om2*erij(k)) + enddo + do k=1,3 + gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)) +!C print *,'gg',k,gg(k) + enddo +! print *,i,j,gg_lipi(3),gg_lipj(3),sss_ele_cut +! write (iout,*) "gg",(gg(k),k=1,3) + do k=1,3 + gvdwx(k,i)=gvdwx(k,i)-gg(k) +gg_lipi(k)& + +(eom12*(dc_norm(k,j)-om12*dc_norm(k,nres+i)) & + +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv + + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)& + +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,j)) & + +eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv + +! write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) & +! +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv +! write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) & +! +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv + enddo +! +! Calculate the components of the gradient in DC and X +! +!grad do k=i,j-1 +!grad do l=1,3 +!grad gvdwc(l,k)=gvdwc(l,k)+gg(l) +!grad enddo +!grad enddo + do l=1,3 + gvdwc(l,i)=gvdwc(l,i)-gg(l) + gvdwc(l,j)=gvdwc(l,j)+gg(l) + enddo + end subroutine sc_grad_cat + + +#ifdef CRYST_THETA +!----------------------------------------------------------------------------- + subroutine mixder(thetai,thet_pred_mean,theta0i,E_tc_t) + + use comm_calcthet +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +! include 'COMMON.LOCAL' +! include 'COMMON.IOUNITS' +!el real(kind=8) :: term1,term2,termm,diffak,ratak,& +!el ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq,& +!el delthe0,sig0inv,sigtc,sigsqtc,delthec, + real(kind=8) :: thetai,thet_pred_mean,theta0i,E_tc_t + real(kind=8) :: t3,t6,t9,t12,t14,t16,t21,t23,t26,t27,t32,t40 +!el integer :: it +!el common /calcthet/ term1,term2,termm,diffak,ratak,& +!el ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq,& +!el delthe0,sig0inv,sigtc,sigsqtc,delthec,it +!el local variables + + delthec=thetai-thet_pred_mean + delthe0=thetai-theta0i +! "Thank you" to MAPLE (probably spared one day of hand-differentiation). + t3 = thetai-thet_pred_mean + t6 = t3**2 + t9 = term1 + t12 = t3*sigcsq + t14 = t12+t6*sigsqtc + t16 = 1.0d0 + t21 = thetai-theta0i + t23 = t21**2 + t26 = term2 + t27 = t21*t26 + t32 = termexp + t40 = t32**2 + E_tc_t = -((sigcsq+2.D0*t3*sigsqtc)*t9-t14*sigcsq*t3*t16*t9 & + -aktc*sig0inv*t27)/t32+(t14*t9+aktc*t26)/t40 & + *(-t12*t9-ak*sig0inv*t27) + return + end subroutine mixder +#endif +!----------------------------------------------------------------------------- +! cartder.F +!----------------------------------------------------------------------------- subroutine cartder !----------------------------------------------------------------------------- ! This subroutine calculates the derivatives of the consecutive virtual @@ -12240,7 +12860,7 @@ ! call intcartderiv ! call checkintcartgrad call zerograd - aincr=2.0D-5 + aincr=1.0D-7 write(iout,*) 'Calling CHECK_ECARTINT.',aincr nf=0 icall=0 @@ -13496,8 +14116,8 @@ rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj))) - sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj))) - sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj))) + sss_ele_cut=sscale_ele(1.0d0/(rij)) + sss_ele_grad=sscagrad_ele(1.0d0/(rij)) sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj))) if (sss_ele_cut.le.0.0) cycle if (sss.lt.1.0d0) then @@ -13553,7 +14173,7 @@ sigder=fac*sigder fac=rij*fac fac=fac+evdwij*(sss_ele_grad/sss_ele_cut& - /sigma(itypi,itypj)*rij-sss_grad/(1.0-sss)*rij & + *rij-sss_grad/(1.0-sss)*rij & /sigmaii(itypi,itypj)) ! fac=0.0d0 ! Calculate the radial part of the gradient @@ -13787,8 +14407,8 @@ rij=dsqrt(rrij) sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj))) sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj))) - sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj))) - sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj))) + sss_ele_cut=sscale_ele(1.0d0/(rij)) + sss_ele_grad=sscagrad_ele(1.0d0/(rij)) if (sss_ele_cut.le.0.0) cycle if (sss.gt.0.0d0) then @@ -13844,7 +14464,7 @@ sigder=fac*sigder fac=rij*fac fac=fac+evdwij*(sss_ele_grad/sss_ele_cut& - /sigma(itypi,itypj)*rij+sss_grad/sss*rij & + *rij+sss_grad/sss*rij & /sigmaii(itypi,itypj)) ! fac=0.0d0 @@ -16123,20 +16743,66 @@ ! ! Calculate the virtual-bond-angle energy. ! - call ebend(ebe,ethetacnstr) -! ! Calculate the SC local energy. ! call vec_and_deriv call esc(escloc) ! + if (wang.gt.0d0) then + if (tor_mode.eq.0) then + call ebend(ebe) + else +!C ebend kcc is Kubo cumulant clustered rigorous attemp to derive the +!C energy function + call ebend_kcc(ebe) + endif + else + ebe=0.0d0 + endif + ethetacnstr=0.0d0 + if (with_theta_constr) call etheta_constr(ethetacnstr) + +! write(iout,*) "in etotal afer ebe",ipot + +! print *,"Processor",myrank," computed UB" +! +! Calculate the SC local energy. +! + call esc(escloc) +!elwrite(iout,*) "in etotal afer esc",ipot +! print *,"Processor",myrank," computed USC" +! +! Calculate the virtual-bond torsional energy. +! +!d print *,'nterm=',nterm +! if (wtor.gt.0) then +! call etor(etors,edihcnstr) +! else +! etors=0 +! edihcnstr=0 +! endif + if (wtor.gt.0.0d0) then + if (tor_mode.eq.0) then + call etor(etors) + else +!C etor kcc is Kubo cumulant clustered rigorous attemp to derive the +!C energy function + call etor_kcc(etors) + endif + else + etors=0.0d0 + endif + edihcnstr=0.0d0 + if (ndih_constr.gt.0) call etor_constr(edihcnstr) + ! Calculate the virtual-bond torsional energy. ! - call etor(etors,edihcnstr) ! ! 6/23/01 Calculate double-torsional energy ! + if ((wtor_d.gt.0.0d0).and.(tor_mode.eq.0)) then call etor_d(etors_d) + endif ! ! 21/5/07 Calculate local sicdechain correlation energy ! @@ -19827,11 +20493,11 @@ integer :: i,j if(nres.lt.100) then - maxconts=nres + maxconts=10*nres elseif(nres.lt.200) then - maxconts=0.8*nres ! Max. number of contacts per residue + maxconts=10*nres ! Max. number of contacts per residue else - maxconts=0.6*nres ! (maxconts=maxres/4) + maxconts=10*nres ! (maxconts=maxres/4) endif maxcont=12*nres ! Max. number of SC contacts maxvar=6*nres ! Max. number of variables @@ -19958,6 +20624,24 @@ allocate(Ug2DtEUgder(2,2,2,nres)) allocate(DtUg2EUgder(2,2,2,nres)) !(2,2,2,maxres) + allocate(b1(2,nres)) !(2,-maxtor:maxtor) + allocate(b2(2,nres)) !(2,-maxtor:maxtor) + allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor) + allocate(b2tilde(2,nres)) !(2,-maxtor:maxtor) + + allocate(ctilde(2,2,nres)) + allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor) + allocate(gtb1(2,nres)) + allocate(gtb2(2,nres)) + allocate(cc(2,2,nres)) + allocate(dd(2,2,nres)) + allocate(ee(2,2,nres)) + allocate(gtcc(2,2,nres)) + allocate(gtdd(2,2,nres)) + allocate(gtee(2,2,nres)) + allocate(gUb2(2,nres)) + allocate(gteUg(2,2,nres)) + ! common /rotat_old/ allocate(costab(nres)) allocate(sintab(nres)) @@ -21188,13 +21872,13 @@ enddo ! IF ( (wcorr_nucl.gt.0.0d0.or.wcorr3_nucl.gt.0.0d0) .and. IF ( j.gt.i+1 .and.& - num_conti.le.maxconts) THEN + num_conti.le.maxcont) THEN !C !C Calculate the contact function. The ith column of the array JCONT will !C contain the numbers of atoms that make contacts with the atom I (of numbers !C greater than I). The arrays FACONT and GACONT will contain the values of !C the contact function and its derivative. - r0ij=2.20D0*sigma(itypi,itypj) + r0ij=2.20D0*sigma_nucl(itypi,itypj) !c write (2,*) "ij",i,j," rij",1.0d0/rij," r0ij",r0ij call gcont(rij,r0ij,1.0D0,0.2d0/r0ij,fcont,fprimcont) !c write (2,*) "fcont",fcont @@ -21204,7 +21888,7 @@ if (num_conti.gt.maxconts) then write (iout,*) 'WARNING - max. # of contacts exceeded;',& - ' will skip next contacts for this conf.' + ' will skip next contacts for this conf.',maxconts else jcont_hb(num_conti,i)=j !c write (iout,*) "num_conti",num_conti, @@ -22067,81 +22751,93 @@ return end subroutine ecatcat !--------------------------------------------------------------------------- - subroutine ecat_prot(ecation_prot) - integer i,j,k,subchap,itmp,inum - real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,& - r7,r4,ecationcation - real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, & - dist_init,dist_temp,ecation_prot,rcal,rocal, & - Evan1,Evan2,EC,cm1mag,DASGL,delta,r0p,Epepcat, & - catl,cml,calpl, Etotal_p, Etotal_m,rtab,wdip,wmodquad,wquad1, & - wquad2,wvan1,E1,E2,wconst,wvan2,rcpm,dcmag,sin2thet,sinthet, & - costhet,v1m,v2m,wh2o,wc,rsecp,Ir,Irsecp,Irthrp,Irfourp,Irfiftp,& - Irsistp,Irseven,Irtwelv,Irthir,dE1dr,dE2dr,dEdcos,wquad2p,opt, & - rs,rthrp,rfourp,rsixp,reight,Irsixp,Ireight,Irtw,Irfourt, & - opt1,opt2,opt3,opt4,opt5,opt6,opt7,opt8,opt9,opt10,opt11,opt12,& - opt13,opt14,opt15,opt16,opt17,opt18,opt19, & - Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip - real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,& - gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, & - dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, & - tab1,tab2,tab3,diff,cm1,sc,p,tcat,talp,cm,drcp,drcp_norm,vcat, & - v1,v2,v3,myd_norm,dx,vcm,valpha,drdpep,dcosdpep,dcosddci,dEdpep,& - dEcCat,dEdipCm,dEdipCalp,dEquad1Cat,dEquad1Cm,dEquad1Calp, & - dEquad2Cat,dEquad2Cm,dEquad2Calpd,Evan1Cat,dEvan1Calp,dEvan2Cat,& - dEvan2Calp,dEtotalCat,dscvec,dEcCm,dEcCalp,dEdipCat,dEquad2Calp,& - dEvan1Cat - real(kind=8),dimension(6) :: vcatprm - ecation_prot=0.0d0 -! first lets calculate interaction with peptide groups - if (nres_molec(5).eq.0) return - 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 +! new for K+ + subroutine ecats_prot_amber(ecations_prot_amber) +! subroutine ecat_prot2(ecation_prot) + use calc_data + use comm_momo + + logical :: lprn +!el local variables + integer :: iint,itypi1,subchap,isel,itmp + real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,e1,e2,sigm,epsi + real(kind=8) :: evdw + real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,& + dist_temp, dist_init,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,egb + real(kind=8) :: Fisocav,ECL,Elj,Equad,Epol,eheadtail,& + Lambf,& + Chif,ChiLambf,Fcav,dFdR,dFdOM1,& + ecations_prot_amber,dFdOM2,dFdL,dFdOM12,& + federmaus,& + d1i,d1j +! real(kind=8),dimension(3,2)::erhead_tail +! real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead,Rtail_distance + real(kind=8) :: facd4, adler, Fgb, facd3 + integer troll,jj,istate + real (kind=8) :: dcosom1(3),dcosom2(3) + + ecations_prot_amber=0.0D0 + if (nres_molec(5).eq.0) return + eps_out=80.0d0 +! sss_ele_cut=1.0d0 + itmp=0 do i=1,4 itmp=itmp+nres_molec(i) enddo ! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization do i=ibond_start,ibond_end -! cycle - if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle ! leave dummy atoms - xi=0.5d0*(c(1,i)+c(1,i+1)) - yi=0.5d0*(c(2,i)+c(2,i+1)) - zi=0.5d0*(c(3,i)+c(3,i+1)) - xi=mod(xi,boxxsize) + +! print *,"I am in EVDW",i + itypi=iabs(itype(i,1)) +! if (i.ne.47) cycle + if (itypi.eq.ntyp1) cycle + itypi1=iabs(itype(i+1,1)) + xi=c(1,nres+i) + yi=c(2,nres+i) + zi=c(3,nres+i) + xi=dmod(xi,boxxsize) if (xi.lt.0) xi=xi+boxxsize - yi=mod(yi,boxysize) + yi=dmod(yi,boxysize) if (yi.lt.0) yi=yi+boxysize - zi=mod(zi,boxzsize) + zi=dmod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize - + dxi=dc_norm(1,nres+i) + dyi=dc_norm(2,nres+i) + dzi=dc_norm(3,nres+i) + dsci_inv=vbld_inv(i+nres) do j=itmp+1,itmp+nres_molec(5) + +! Calculate SC interaction energy. + itypj=iabs(itype(j,1)) + if ((itypj.eq.ntyp1)) cycle + CALL elgrad_init_cat(eheadtail,Egb,Ecl,Elj,Equad,Epol) + + dscj_inv=vbld_inv(j+nres) xj=c(1,j) yj=c(2,j) zj=c(3,j) - xj=dmod(xj,boxxsize) - if (xj.lt.0) xj=xj+boxxsize - yj=dmod(yj,boxysize) - if (yj.lt.0) yj=yj+boxysize - zj=dmod(zj,boxzsize) - if (zj.lt.0) zj=zj+boxzsize - dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 - xj_safe=xj - yj_safe=yj - zj_safe=zj - subchap=0 - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 + xj=dmod(xj,boxxsize) + if (xj.lt.0) xj=xj+boxxsize + yj=dmod(yj,boxysize) + if (yj.lt.0) yj=yj+boxysize + zj=dmod(zj,boxzsize) + if (zj.lt.0) zj=zj+boxzsize + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 + + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 xj=xj_safe+xshift*boxxsize yj=yj_safe+yshift*boxysize zj=zj_safe+zshift*boxzsize @@ -22153,30 +22849,422 @@ zj_temp=zj subchap=1 endif - enddo - enddo - enddo - if (subchap.eq.1) then + enddo + enddo + enddo + if (subchap.eq.1) then xj=xj_temp-xi yj=yj_temp-yi zj=zj_temp-zi - else + else xj=xj_safe-xi yj=yj_safe-yi zj=zj_safe-zi - endif -! enddo -! enddo - rcpm = sqrt(xj**2+yj**2+zj**2) - drcp_norm(1)=xj/rcpm - drcp_norm(2)=yj/rcpm - drcp_norm(3)=zj/rcpm - dcmag=0.0 - do k=1,3 - dcmag=dcmag+dc(k,i)**2 - enddo - dcmag=dsqrt(dcmag) - do k=1,3 + endif + +! dxj = dc_norm( 1, nres+j ) +! dyj = dc_norm( 2, nres+j ) +! dzj = dc_norm( 3, nres+j ) + + itypi = itype(i,1) + itypj = itype(j,5) +! Parameters from fitting the analitical expressions to the PMF obtained by umbrella +! sampling performed with amber package +! alf1 = 0.0d0 +! alf2 = 0.0d0 +! alf12 = 0.0d0 +! a12sq = rborn(itypi,itypj) * rborn(itypj,itypi) + chi1 = chicat(itypi,itypj) + chis1 = chiscat(itypi,itypj) + chip1 = chippcat(itypi,itypj) +! chis2 = chis(itypj,itypi) +! chis12 = chis1 * chis2 + sig1 = sigmap1cat(itypi,itypj) +! sig2 = sigmap2(itypi,itypj) +! alpha factors from Fcav/Gcav + b1cav = alphasurcat(1,itypi,itypj) + b2cav = alphasurcat(2,itypi,itypj) + b3cav = alphasurcat(3,itypi,itypj) + b4cav = alphasurcat(4,itypi,itypj) + +! used to determine whether we want to do quadrupole calculations + eps_in = epsintabcat(itypi,itypj) + if (eps_in.eq.0.0) eps_in=1.0 + + eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out)) +! Rtail = 0.0d0 + + DO k = 1, 3 + ctail(k,1)=c(k,i+nres) + ctail(k,2)=c(k,j) + END DO +!c! tail distances will be themselves usefull elswhere +!c1 (in Gcav, for example) + Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 ) + Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 ) + Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 ) + Rtail = dsqrt( & + (Rtail_distance(1)*Rtail_distance(1)) & + + (Rtail_distance(2)*Rtail_distance(2)) & + + (Rtail_distance(3)*Rtail_distance(3))) +! tail location and distance calculations +! dhead1 + d1 = dheadcat(1, 1, itypi, itypj) +! d2 = dhead(2, 1, itypi, itypj) + DO k = 1,3 +! location of polar head is computed by taking hydrophobic centre +! and moving by a d1 * dc_norm vector +! see unres publications for very informative images + chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres) + chead(k,2) = c(k, j) +! distance +! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres)) +! Rsc(k) = Rsc_distance(k) * Rsc_distance(k) + Rhead_distance(k) = chead(k,2) - chead(k,1) + END DO +! pitagoras (root of sum of squares) + Rhead = dsqrt( & + (Rhead_distance(1)*Rhead_distance(1)) & + + (Rhead_distance(2)*Rhead_distance(2)) & + + (Rhead_distance(3)*Rhead_distance(3))) +!------------------------------------------------------------------- +! zero everything that should be zero'ed + evdwij = 0.0d0 + ECL = 0.0d0 + Elj = 0.0d0 + Equad = 0.0d0 + Epol = 0.0d0 + Fcav=0.0d0 + eheadtail = 0.0d0 + dGCLdOM1 = 0.0d0 + dGCLdOM2 = 0.0d0 + dGCLdOM12 = 0.0d0 + dPOLdOM1 = 0.0d0 + dPOLdOM2 = 0.0d0 + Fcav = 0.0d0 + dFdR = 0.0d0 + dCAVdOM1 = 0.0d0 + dCAVdOM2 = 0.0d0 + dCAVdOM12 = 0.0d0 + dscj_inv = vbld_inv(j+nres) +! print *,i,j,dscj_inv,dsci_inv +! rij holds 1/(distance of Calpha atoms) + rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj) + rij = dsqrt(rrij) + CALL sc_angular +! this should be in elgrad_init but om's are calculated by sc_angular +! which in turn is used by older potentials +! om = omega, sqom = om^2 + sqom1 = om1 * om1 + sqom2 = om2 * om2 + sqom12 = om12 * om12 + +! now we calculate EGB - Gey-Berne +! It will be summed up in evdwij and saved in evdw + sigsq = 1.0D0 / sigsq + sig = sig0ij * dsqrt(sigsq) +! rij_shift = 1.0D0 / rij - sig + sig0ij + rij_shift = Rtail - sig + sig0ij + IF (rij_shift.le.0.0D0) THEN + evdw = 1.0D20 + RETURN + END IF + sigder = -sig * sigsq + rij_shift = 1.0D0 / rij_shift + fac = rij_shift**expon + c1 = fac * fac * aa_aq(itypi,itypj) +! print *,"ADAM",aa_aq(itypi,itypj) + +! c1 = 0.0d0 + c2 = fac * bb_aq(itypi,itypj) +! c2 = 0.0d0 + evdwij = eps1 * eps2rt * eps3rt * ( c1 + c2 ) + eps2der = eps3rt * evdwij + eps3der = eps2rt * evdwij +! evdwij = 4.0d0 * eps2rt * eps3rt * evdwij + evdwij = eps2rt * eps3rt * evdwij +!#ifdef TSCSC +! IF (bb_aq(itypi,itypj).gt.0) THEN +! evdw_p = evdw_p + evdwij +! ELSE +! evdw_m = evdw_m + evdwij +! END IF +!#else + evdw = evdw & + + evdwij +!#endif + c1 = c1 * eps1 * eps2rt**2 * eps3rt**2 + fac = -expon * (c1 + evdwij) * rij_shift + sigder = fac * sigder +! Calculate distance derivative + gg(1) = fac + gg(2) = fac + gg(3) = fac + + fac = chis1 * sqom1 + chis2 * sqom2 & + - 2.0d0 * chis12 * om1 * om2 * om12 + pom = 1.0d0 - chis1 * chis2 * sqom12 + Lambf = (1.0d0 - (fac / pom)) + Lambf = dsqrt(Lambf) + sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0) + Chif = Rtail * sparrow + ChiLambf = Chif * Lambf + eagle = dsqrt(ChiLambf) + bat = ChiLambf ** 11.0d0 + top = b1cav * ( eagle + b2cav * ChiLambf - b3cav ) + bot = 1.0d0 + b4cav * (ChiLambf ** 12.0d0) + botsq = bot * bot + Fcav = top / bot + + dtop = b1cav * ((Lambf / (2.0d0 * eagle)) + (b2cav * Lambf)) + dbot = 12.0d0 * b4cav * bat * Lambf + dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow + + dtop = b1cav * ((Chif / (2.0d0 * eagle)) + (b2cav * Chif)) + dbot = 12.0d0 * b4cav * bat * Chif + eagle = Lambf * pom + dFdOM1 = -(chis1 * om1 - chis12 * om2 * om12) / (eagle) + dFdOM2 = -(chis2 * om2 - chis12 * om1 * om12) / (eagle) + dFdOM12 = chis12 * (chis1 * om1 * om12 - om2) & + * (chis2 * om2 * om12 - om1) / (eagle * pom) + + dFdL = ((dtop * bot - top * dbot) / botsq) + dCAVdOM1 = dFdL * ( dFdOM1 ) + dCAVdOM2 = dFdL * ( dFdOM2 ) + dCAVdOM12 = dFdL * ( dFdOM12 ) + + DO k= 1, 3 + ertail(k) = Rtail_distance(k)/Rtail + END DO + erdxi = scalar( ertail(1), dC_norm(1,i+nres) ) + erdxj = scalar( ertail(1), dC_norm(1,j) ) + facd1 = dtail(1,itypi,itypj) * vbld_inv(i+nres) + facd2 = dtail(2,itypi,itypj) * vbld_inv(j+nres) + DO k = 1, 3 + pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres)) + gvdwx(k,i) = gvdwx(k,i) & + - (( dFdR + gg(k) ) * pom) + pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres)) +! gvdwx(k,j) = gvdwx(k,j) & +! + (( dFdR + gg(k) ) * pom) + gvdwc(k,i) = gvdwc(k,i) & + - (( dFdR + gg(k) ) * ertail(k)) + gvdwc(k,j) = gvdwc(k,j) & + + (( dFdR + gg(k) ) * ertail(k)) + gg(k) = 0.0d0 + +!c! Compute head-head and head-tail energies for each state + isel = iabs(Qi) + iabs(Qj) + IF (isel.eq.0) THEN +!c! No charges - do nothing + eheadtail = 0.0d0 + + ELSE IF (isel.eq.1 .and. iabs(Qj).eq.1) THEN +!c! Nonpolar-charge interactions + if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then + Qi=Qi*2 + Qij=Qij*2 + endif + if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then + Qj=Qj*2 + Qij=Qij*2 + endif + + CALL enq_cat(epol) + eheadtail = epol +! eheadtail = 0.0d0 + + ELSE IF (isel.eq.3 .and. icharge(itypj).eq.2) THEN +!c! Dipole-charge interactions + if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then + Qi=Qi*2 + Qij=Qij*2 + endif + if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then + Qj=Qj*2 + Qij=Qij*2 + endif + CALL edq_cat(ecl, elj, epol) + eheadtail = ECL + elj + epol +! eheadtail = 0.0d0 + + ELSE IF ((isel.eq.2.and. & + iabs(Qi).eq.1).and. & + nstate(itypi,itypj).eq.1) THEN + +!c! Same charge-charge interaction ( +/+ or -/- ) + if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then + Qi=Qi*2 + Qij=Qij*2 + endif + if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then + Qj=Qj*2 + Qij=Qij*2 + endif + + CALL eqq_cat(Ecl,Egb,Epol,Fisocav,Elj) + eheadtail = ECL + Egb + Epol + Fisocav + Elj +! eheadtail = 0.0d0 + +! ELSE IF ((isel.eq.2.and. & +! iabs(Qi).eq.1).and. & +! nstate(itypi,itypj).ne.1) THEN +!c! Different charge-charge interaction ( +/- or -/+ ) +! if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then +! Qi=Qi*2 +! Qij=Qij*2 +! endif +! if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then +! Qj=Qj*2 +! Qij=Qij*2 +! endif +! +! CALL energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad) + END IF ! this endif ends the "catch the gly-gly" at the beggining of Fcav + evdw = evdw + Fcav + eheadtail + + IF (energy_dec) write (iout,'(2(1x,a3,i3),3f6.2,10f16.7)') & + restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& + 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& + Equad,evdwij+Fcav+eheadtail,evdw +! evdw = evdw + Fcav + eheadtail + +! iF (nstate(itypi,itypj).eq.1) THEN + CALL sc_grad_cat +! END IF +!c!------------------------------------------------------------------- +!c! NAPISY KONCOWE + END DO ! j + END DO ! iint + END DO ! i +!c write (iout,*) "Number of loop steps in EGB:",ind +!c energy_dec=.false. +! print *,"EVDW KURW",evdw,nres + + return + end subroutine ecats_prot_amber + +!--------------------------------------------------------------------------- +! old for Ca2+ + subroutine ecat_prot(ecation_prot) +! use calc_data +! use comm_momo + integer i,j,k,subchap,itmp,inum + real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,& + r7,r4,ecationcation + real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, & + dist_init,dist_temp,ecation_prot,rcal,rocal, & + Evan1,Evan2,EC,cm1mag,DASGL,delta,r0p,Epepcat, & + catl,cml,calpl, Etotal_p, Etotal_m,rtab,wdip,wmodquad,wquad1, & + wquad2,wvan1,E1,E2,wconst,wvan2,rcpm,dcmag,sin2thet,sinthet, & + costhet,v1m,v2m,wh2o,wc,rsecp,Ir,Irsecp,Irthrp,Irfourp,Irfiftp,& + Irsistp,Irseven,Irtwelv,Irthir,dE1dr,dE2dr,dEdcos,wquad2p,opt, & + rs,rthrp,rfourp,rsixp,reight,Irsixp,Ireight,Irtw,Irfourt, & + opt1,opt2,opt3,opt4,opt5,opt6,opt7,opt8,opt9,opt10,opt11,opt12,& + opt13,opt14,opt15,opt16,opt17,opt18,opt19, & + Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip,& + ndiv,ndivi + real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,& + gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, & + dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, & + tab1,tab2,tab3,diff,cm1,sc,p,tcat,talp,cm,drcp,drcp_norm,vcat, & + v1,v2,v3,myd_norm,dx,vcm,valpha,drdpep,dcosdpep,dcosddci,dEdpep,& + dEcCat,dEdipCm,dEdipCalp,dEquad1Cat,dEquad1Cm,dEquad1Calp, & + dEquad2Cat,dEquad2Cm,dEquad2Calpd,Evan1Cat,dEvan1Calp,dEvan2Cat,& + dEvan2Calp,dEtotalCat,dscvec,dEcCm,dEcCalp,dEdipCat,dEquad2Calp,& + dEvan1Cat + real(kind=8),dimension(6) :: vcatprm + ecation_prot=0.0d0 +! first lets calculate interaction with peptide groups + if (nres_molec(5).eq.0) return + itmp=0 + do i=1,4 + itmp=itmp+nres_molec(i) + enddo +! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization + do i=ibond_start,ibond_end +! cycle + if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle ! leave dummy atoms + xi=0.5d0*(c(1,i)+c(1,i+1)) + yi=0.5d0*(c(2,i)+c(2,i+1)) + zi=0.5d0*(c(3,i)+c(3,i+1)) + xi=mod(xi,boxxsize) + if (xi.lt.0) xi=xi+boxxsize + yi=mod(yi,boxysize) + if (yi.lt.0) yi=yi+boxysize + zi=mod(zi,boxzsize) + if (zi.lt.0) zi=zi+boxzsize + + do j=itmp+1,itmp+nres_molec(5) +! print *,"WTF",itmp,j,i +! all parameters were for Ca2+ to approximate single charge divide by two + ndiv=1.0 + if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0 + wconst=78*ndiv + wdip =1.092777950857032D2 + wdip=wdip/wconst + wmodquad=-2.174122713004870D4 + wmodquad=wmodquad/wconst + wquad1 = 3.901232068562804D1 + wquad1=wquad1/wconst + wquad2 = 3 + wquad2=wquad2/wconst + wvan1 = 0.1 + wvan2 = 6 +! itmp=0 + + xj=c(1,j) + yj=c(2,j) + zj=c(3,j) + xj=dmod(xj,boxxsize) + if (xj.lt.0) xj=xj+boxxsize + yj=dmod(yj,boxysize) + if (yj.lt.0) yj=yj+boxysize + zj=dmod(zj,boxzsize) + if (zj.lt.0) zj=zj+boxzsize + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif +! enddo +! enddo + rcpm = sqrt(xj**2+yj**2+zj**2) + drcp_norm(1)=xj/rcpm + drcp_norm(2)=yj/rcpm + drcp_norm(3)=zj/rcpm + dcmag=0.0 + do k=1,3 + dcmag=dcmag+dc(k,i)**2 + enddo + dcmag=dsqrt(dcmag) + do k=1,3 myd_norm(k)=dc(k,i)/dcmag enddo costhet=drcp_norm(1)*myd_norm(1)+drcp_norm(2)*myd_norm(2)+& @@ -22198,6 +23286,7 @@ E2 = -wquad1*Irthrp*wquad2+wvan1*(wvan2**12*Irtwelv-& 2*wvan2**6*Irsistp) ecation_prot = ecation_prot+E1+E2 +! print *,"ecatprot",i,j,ecation_prot,rcpm dE1dr = -2*costhet*wdip*Irthrp-& (4*wmodquad*Irfiftp+3*wquad1*Irfourp)*sin2thet dE2dr = 3*wquad1*wquad2*Irfourp- & @@ -22237,6 +23326,9 @@ enddo cm1mag=sqrt(cm1(1)**2+cm1(2)**2+cm1(3)**2) do j=itmp+1,itmp+nres_molec(5) + ndiv=1.0 + if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0 + xj=c(1,j) yj=c(2,j) zj=c(3,j) @@ -22279,7 +23371,10 @@ endif ! enddo ! enddo - if(itype(i,1).eq.15.or.itype(i,1).eq.16) then +! 15- Glu 16-Asp + if((itype(i,1).eq.15.or.itype(i,1).eq.16).or.& + ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.& + (itype(i,1).eq.25))) then if(itype(i,1).eq.16) then inum=1 else @@ -22289,12 +23384,29 @@ 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 +! do k=1,3 +! vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres) + vcm(1)=(cm1(1)/cm1mag)*dASGL+xi + vcm(2)=(cm1(2)/cm1mag)*dASGL+yi + vcm(3)=(cm1(3)/cm1mag)*dASGL+zi + +! valpha(k)=c(k,i) +! vcat(k)=c(k,j) + if (subchap.eq.1) then + vcat(1)=xj_temp + vcat(2)=yj_temp + vcat(3)=zj_temp + else + vcat(1)=xj_safe + vcat(2)=yj_safe + vcat(3)=zj_safe + endif + valpha(1)=xi-c(1,i+nres)+c(1,i) + valpha(2)=yi-c(2,i+nres)+c(2,i) + valpha(3)=zi-c(3,i+nres)+c(3,i) + +! enddo + do k=1,3 dx(k) = vcat(k)-vcm(k) enddo do k=1,3 @@ -22307,7 +23419,15 @@ ! The weights of the energy function calculated from !The quantum mechanical GAMESS simulations of calcium with ASP/GLU - wh2o=78 + if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then + ndivi=0.5 + else + ndivi=1.0 + endif + ndiv=1.0 + if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0 + + wh2o=78*ndivi*ndiv wc = vcatprm(1) wc=wc/wh2o wdip =vcatprm(2) @@ -22316,7 +23436,7 @@ wquad1=wquad1/wh2o wquad2 = vcatprm(4) wquad2=wquad2/wh2o - wquad2p = 1-wquad2 + wquad2p = 1.0d0-wquad2 wvan1 = vcatprm(5) wvan2 =vcatprm(6) opt = dx(1)**2+dx(2)**2 @@ -22327,11 +23447,11 @@ rsixp = rfourp*rsecp reight=rsixp*rsecp Ir = 1.0d0/rs - Irsecp = 1/rsecp + Irsecp = 1.0d0/rsecp Irthrp = Irsecp/rs Irfourp = Irthrp/rs - Irsixp = 1/rsixp - Ireight=1/reight + Irsixp = 1.0d0/rsixp + Ireight=1.0d0/reight Irtw=Irsixp*Irsixp Irthir=Irtw/rs Irfourt=Irthir/rs @@ -22442,11 +23562,27 @@ 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 +! vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres) +! valpha(k)=c(k,i) +! vcat(k)=c(k,j) +! enddo + vcm(1)=(cm1(1)/cm1mag)*dASGL+xi + vcm(2)=(cm1(2)/cm1mag)*dASGL+yi + vcm(3)=(cm1(3)/cm1mag)*dASGL+zi + if (subchap.eq.1) then + vcat(1)=xj_temp + vcat(2)=yj_temp + vcat(3)=zj_temp + else + vcat(1)=xj_safe + vcat(2)=yj_safe + vcat(3)=zj_safe + endif + valpha(1)=xi-c(1,i+nres)+c(1,i) + valpha(2)=yi-c(2,i+nres)+c(2,i) + valpha(3)=zi-c(3,i+nres)+c(3,i) + do k=1,3 dx(k) = vcat(k)-vcm(k) @@ -22460,7 +23596,10 @@ 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 + ndiv=1.0 + if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0 + + wh2o=78*ndiv wdip =vcatprm(2) wdip=wdip/wh2o wquad1 =vcatprm(3) @@ -22555,6 +23694,11 @@ dscmag = 0.0d0 do k=1,3 dscvec(k) = c(k,i+nres)-c(k,i) +! TU SPRAWDZ??? +! dscvec(1) = xj +! dscvec(2) = yj +! dscvec(3) = zj + dscmag = dscmag+dscvec(k)*dscvec(k) enddo dscmag3 = dscmag @@ -22576,7 +23720,10 @@ else rcal = 0.0d0 do k=1,3 - r(k) = c(k,j)-c(k,i+nres) +! r(k) = c(k,j)-c(k,i+nres) + r(1) = xj + r(2) = yj + r(3) = zj rcal = rcal+r(k)*r(k) enddo ract=sqrt(rcal) @@ -24720,6 +25867,16 @@ !c! Compute head-head and head-tail energies for each state isel = iabs(Qi) + iabs(Qj) +! double charge for Phophorylated! itype - 25,27,27 +! if ((itype(i).eq.27).or.(itype(i).eq.26).or.(itype(i).eq.25)) then +! Qi=Qi*2 +! Qij=Qij*2 +! endif +! if ((itype(j).eq.27).or.(itype(j).eq.26).or.(itype(j).eq.25)) then +! Qj=Qj*2 +! Qij=Qij*2 +! endif + ! isel=0 IF (isel.eq.0) THEN !c! No charges - do nothing @@ -24733,24 +25890,59 @@ ELSE IF (isel.eq.1 .and. iabs(Qi).eq.1) THEN !c! Charge-nonpolar interactions + if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then + Qi=Qi*2 + Qij=Qij*2 + endif + if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then + Qj=Qj*2 + Qij=Qij*2 + endif + CALL eqn(epol) eheadtail = epol ! eheadtail = 0.0d0 ELSE IF (isel.eq.1 .and. iabs(Qj).eq.1) THEN !c! Nonpolar-charge interactions + if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then + Qi=Qi*2 + Qij=Qij*2 + endif + if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then + Qj=Qj*2 + Qij=Qij*2 + endif + CALL enq(epol) eheadtail = epol ! eheadtail = 0.0d0 ELSE IF (isel.eq.3 .and. icharge(itypj).eq.2) THEN !c! Charge-dipole interactions + if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then + Qi=Qi*2 + Qij=Qij*2 + endif + if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then + Qj=Qj*2 + Qij=Qij*2 + endif + CALL eqd(ecl, elj, epol) eheadtail = ECL + elj + epol ! eheadtail = 0.0d0 ELSE IF (isel.eq.3 .and. icharge(itypi).eq.2) THEN !c! Dipole-charge interactions + if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then + Qi=Qi*2 + Qij=Qij*2 + endif + if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then + Qj=Qj*2 + Qij=Qij*2 + endif CALL edq(ecl, elj, epol) eheadtail = ECL + elj + epol ! eheadtail = 0.0d0 @@ -24759,6 +25951,15 @@ iabs(Qi).eq.1).and. & nstate(itypi,itypj).eq.1) THEN !c! Same charge-charge interaction ( +/+ or -/- ) + if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then + Qi=Qi*2 + Qij=Qij*2 + endif + if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then + Qj=Qj*2 + Qij=Qij*2 + endif + CALL eqq(Ecl,Egb,Epol,Fisocav,Elj) eheadtail = ECL + Egb + Epol + Fisocav + Elj ! eheadtail = 0.0d0 @@ -24767,6 +25968,15 @@ iabs(Qi).eq.1).and. & nstate(itypi,itypj).ne.1) THEN !c! Different charge-charge interaction ( +/- or -/+ ) + if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then + Qi=Qi*2 + Qij=Qij*2 + endif + if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then + Qj=Qj*2 + Qij=Qij*2 + endif + CALL energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad) END IF END IF ! this endif ends the "catch the gly-gly" at the beggining of Fcav @@ -24797,7 +26007,7 @@ use calc_data use comm_momo real (kind=8) :: facd3, facd4, federmaus, adler,& - Ecl,Egb,Epol,Fisocav,Elj,Fgb + Ecl,Egb,Epol,Fisocav,Elj,Fgb,debkap ! integer :: k !c! Epol and Gpol analytical parameters alphapol1 = alphapol(itypi,itypj) @@ -24842,10 +26052,15 @@ dGCLdOM12 = 0.0d0 ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq)) Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0) - Egb = -(332.0d0 * Qij * eps_inout_fac) / Fgb + debkap=debaykap(itypi,itypj) + Egb = -(332.0d0 * Qij *& + (1.0/eps_in-dexp(-debkap*Fgb)/eps_out)) / Fgb ! print *,"EGB WTF",Qij,eps_inout_fac,Fgb,itypi,itypj,eps_in,eps_out !c! Derivative of Egb is Ggb... - dGGBdFGB = -(-332.0d0 * Qij * eps_inout_fac) / (Fgb * Fgb) + dGGBdFGB = -(-332.0d0 * Qij * & + (1.0/eps_in-dexp(-debkap*Fgb)/eps_out))/(Fgb*Fgb)& + -(332.0d0 * Qij *& + (dexp(-debkap*Fgb)*debkap/eps_out))/ Fgb dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )/ ( 2.0d0 * Fgb ) dGGBdR = dGGBdFGB * dFGBdR !c!------------------------------------------------------------------- @@ -24977,6 +26192,196 @@ END DO RETURN END SUBROUTINE eqq + + SUBROUTINE eqq_cat(Ecl,Egb,Epol,Fisocav,Elj) + use calc_data + use comm_momo + real (kind=8) :: facd3, facd4, federmaus, adler,& + Ecl,Egb,Epol,Fisocav,Elj,Fgb,debkap +! integer :: k +!c! Epol and Gpol analytical parameters + alphapol1 = alphapolcat(itypi,itypj) + alphapol2 = alphapolcat(itypj,itypi) +!c! Fisocav and Gisocav analytical parameters + al1 = alphisocat(1,itypi,itypj) + al2 = alphisocat(2,itypi,itypj) + al3 = alphisocat(3,itypi,itypj) + al4 = alphisocat(4,itypi,itypj) + csig = (1.0d0 & + / dsqrt(sigiso1cat(itypi, itypj)**2.0d0 & + + sigiso2cat(itypi,itypj)**2.0d0)) +!c! + pis = sig0headcat(itypi,itypj) + eps_head = epsheadcat(itypi,itypj) + Rhead_sq = Rhead * Rhead +!c! R1 - distance between head of ith side chain and tail of jth sidechain +!c! R2 - distance between head of jth side chain and tail of ith sidechain + R1 = 0.0d0 + R2 = 0.0d0 + DO k = 1, 3 +!c! Calculate head-to-tail distances needed by Epol + R1=R1+(ctail(k,2)-chead(k,1))**2 + R2=R2+(chead(k,2)-ctail(k,1))**2 + END DO +!c! Pitagoras + R1 = dsqrt(R1) + R2 = dsqrt(R2) + +!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj) +!c! & +dhead(1,1,itypi,itypj))**2)) +!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj) +!c! & +dhead(2,1,itypi,itypj))**2)) + +!c!------------------------------------------------------------------- +!c! Coulomb electrostatic interaction + Ecl = (332.0d0 * Qij) / Rhead +!c! derivative of Ecl is Gcl... + dGCLdR = (-332.0d0 * Qij ) / Rhead_sq + dGCLdOM1 = 0.0d0 + dGCLdOM2 = 0.0d0 + dGCLdOM12 = 0.0d0 + ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq)) + Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0) + debkap=debaykapcat(itypi,itypj) + Egb = -(332.0d0 * Qij *& + (1.0/eps_in-dexp(-debkap*Fgb)/eps_out)) / Fgb +! print *,"EGB WTF",Qij,eps_inout_fac,Fgb,itypi,itypj,eps_in,eps_out +!c! Derivative of Egb is Ggb... + dGGBdFGB = -(-332.0d0 * Qij * & + (1.0/eps_in-dexp(-debkap*Fgb)/eps_out))/(Fgb*Fgb)& + -(332.0d0 * Qij *& + (dexp(-debkap*Fgb)*debkap/eps_out))/ Fgb + dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )/ ( 2.0d0 * Fgb ) + dGGBdR = dGGBdFGB * dFGBdR +!c!------------------------------------------------------------------- +!c! Fisocav - isotropic cavity creation term +!c! or "how much energy it costs to put charged head in water" + pom = Rhead * csig + top = al1 * (dsqrt(pom) + al2 * pom - al3) + bot = (1.0d0 + al4 * pom**12.0d0) + botsq = bot * bot + FisoCav = top / bot +! write (*,*) "Rhead = ",Rhead +! write (*,*) "csig = ",csig +! write (*,*) "pom = ",pom +! write (*,*) "al1 = ",al1 +! write (*,*) "al2 = ",al2 +! write (*,*) "al3 = ",al3 +! write (*,*) "al4 = ",al4 +! write (*,*) "top = ",top +! write (*,*) "bot = ",bot +!c! Derivative of Fisocav is GCV... + dtop = al1 * ((1.0d0 / (2.0d0 * dsqrt(pom))) + al2) + dbot = 12.0d0 * al4 * pom ** 11.0d0 + dGCVdR = ((dtop * bot - top * dbot) / botsq) * csig +!c!------------------------------------------------------------------- +!c! Epol +!c! Polarization energy - charged heads polarize hydrophobic "neck" + MomoFac1 = (1.0d0 - chi1 * sqom2) + MomoFac2 = (1.0d0 - chi2 * sqom1) + RR1 = ( R1 * R1 ) / MomoFac1 + RR2 = ( R2 * R2 ) / MomoFac2 + ee1 = exp(-( RR1 / (4.0d0 * a12sq) )) + ee2 = exp(-( RR2 / (4.0d0 * a12sq) )) + fgb1 = sqrt( RR1 + a12sq * ee1 ) + fgb2 = sqrt( RR2 + a12sq * ee2 ) + epol = 332.0d0 * eps_inout_fac * ( & + (( alphapol1 / fgb1 )**4.0d0)+((alphapol2/fgb2) ** 4.0d0 )) +!c! epol = 0.0d0 + dPOLdFGB1 = -(1328.0d0 * eps_inout_fac * alphapol1 ** 4.0d0)& + / (fgb1 ** 5.0d0) + dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0)& + / (fgb2 ** 5.0d0) + dFGBdR1 = ( (R1 / MomoFac1)* ( 2.0d0 - (0.5d0 * ee1) ) )& + / ( 2.0d0 * fgb1 ) + dFGBdR2 = ( (R2 / MomoFac2)* ( 2.0d0 - (0.5d0 * ee2) ) )& + / ( 2.0d0 * fgb2 ) + dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1))& + * ( 2.0d0 - 0.5d0 * ee1) ) / ( 2.0d0 * fgb1 ) + dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2))& + * ( 2.0d0 - 0.5d0 * ee2) ) / ( 2.0d0 * fgb2 ) + dPOLdR1 = dPOLdFGB1 * dFGBdR1 +!c! dPOLdR1 = 0.0d0 + dPOLdR2 = dPOLdFGB2 * dFGBdR2 +!c! dPOLdR2 = 0.0d0 + dPOLdOM1 = dPOLdFGB2 * dFGBdOM1 +!c! dPOLdOM1 = 0.0d0 + dPOLdOM2 = dPOLdFGB1 * dFGBdOM2 +!c! dPOLdOM2 = 0.0d0 +!c!------------------------------------------------------------------- +!c! Elj +!c! Lennard-Jones 6-12 interaction between heads + pom = (pis / Rhead)**6.0d0 + Elj = 4.0d0 * eps_head * pom * (pom-1.0d0) +!c! derivative of Elj is Glj + dGLJdR = 4.0d0 * eps_head*(((-12.0d0*pis**12.0d0)/(Rhead**13.0d0))& + + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0))) +!c!------------------------------------------------------------------- +!c! Return the results +!c! These things do the dRdX derivatives, that is +!c! allow us to change what we see from function that changes with +!c! distance to function that changes with LOCATION (of the interaction +!c! site) + DO k = 1, 3 + erhead(k) = Rhead_distance(k)/Rhead + erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1) + erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2) + END DO + + erdxi = scalar( erhead(1), dC_norm(1,i+nres) ) + erdxj = scalar( erhead(1), dC_norm(1,j) ) + bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) ) + federmaus = scalar(erhead_tail(1,1),dC_norm(1,j)) + eagle = scalar( erhead_tail(1,2), dC_norm(1,j) ) + adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) ) + facd1 = d1 * vbld_inv(i+nres) + facd2 = d2 * vbld_inv(j) + facd3 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres) + facd4 = dtailcat(2,itypi,itypj) * vbld_inv(j) + +!c! Now we add appropriate partial derivatives (one in each dimension) + DO k = 1, 3 + hawk = (erhead_tail(k,1) + & + facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres))) + condor = (erhead_tail(k,2) + & + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j))) + + pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres)) + gvdwx(k,i) = gvdwx(k,i) & + - dGCLdR * pom& + - dGGBdR * pom& + - dGCVdR * pom& + - dPOLdR1 * hawk& + - dPOLdR2 * (erhead_tail(k,2)& + -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))& + - dGLJdR * pom + + pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j)) + gvdwx(k,j) = gvdwx(k,j)+ dGCLdR * pom& + + dGGBdR * pom+ dGCVdR * pom& + + dPOLdR1 * (erhead_tail(k,1)& + -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j)))& + + dPOLdR2 * condor + dGLJdR * pom + + gvdwc(k,i) = gvdwc(k,i) & + - dGCLdR * erhead(k)& + - dGGBdR * erhead(k)& + - dGCVdR * erhead(k)& + - dPOLdR1 * erhead_tail(k,1)& + - dPOLdR2 * erhead_tail(k,2)& + - dGLJdR * erhead(k) + + gvdwc(k,j) = gvdwc(k,j) & + + dGCLdR * erhead(k) & + + dGGBdR * erhead(k) & + + dGCVdR * erhead(k) & + + dPOLdR1 * erhead_tail(k,1) & + + dPOLdR2 * erhead_tail(k,2)& + + dGLJdR * erhead(k) + + END DO + RETURN + END SUBROUTINE eqq_cat !c!------------------------------------------------------------------- SUBROUTINE energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad) use comm_momo @@ -25363,26 +26768,93 @@ 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))) + hawk = (erhead_tail(k,1) + & + facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres))) + + gvdwx(k,i) = gvdwx(k,i) & + - dPOLdR1 * hawk + gvdwx(k,j) = gvdwx(k,j) & + + dPOLdR1 * (erhead_tail(k,1) & + -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) + + gvdwc(k,i) = gvdwc(k,i) - dPOLdR1 * erhead_tail(k,1) + gvdwc(k,j) = gvdwc(k,j) + dPOLdR1 * erhead_tail(k,1) + + END DO + RETURN + END SUBROUTINE eqn + SUBROUTINE enq(Epol) + use calc_data + use comm_momo + double precision facd3, adler,epol + alphapol2 = alphapol(itypj,itypi) +!c! R2 - distance between head of jth side chain and tail of ith sidechain + R2 = 0.0d0 + DO k = 1, 3 +!c! Calculate head-to-tail distances + R2=R2+(chead(k,2)-ctail(k,1))**2 + END DO +!c! Pitagoras + R2 = dsqrt(R2) + +!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj) +!c! & +dhead(1,1,itypi,itypj))**2)) +!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj) +!c! & +dhead(2,1,itypi,itypj))**2)) +!c------------------------------------------------------------------------ +!c Polarization energy + MomoFac2 = (1.0d0 - chi2 * sqom1) + RR2 = R2 * R2 / MomoFac2 + ee2 = exp(-(RR2 / (4.0d0 * a12sq))) + fgb2 = sqrt(RR2 + a12sq * ee2) + epol = 332.0d0 * eps_inout_fac * ((alphapol2/fgb2) ** 4.0d0 ) + dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0) & + / (fgb2 ** 5.0d0) + dFGBdR2 = ( (R2 / MomoFac2) & + * ( 2.0d0 - (0.5d0 * ee2) ) ) & + / (2.0d0 * fgb2) + dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) & + * (2.0d0 - 0.5d0 * ee2) ) & + / (2.0d0 * fgb2) + dPOLdR2 = dPOLdFGB2 * dFGBdR2 +!c! dPOLdR2 = 0.0d0 + dPOLdOM1 = dPOLdFGB2 * dFGBdOM1 +!c! dPOLdOM1 = 0.0d0 + dPOLdOM2 = 0.0d0 +!c!------------------------------------------------------------------- +!c! Return the results +!c! (See comments in Eqq) + DO k = 1, 3 + erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2) + END DO + eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) ) + adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) ) + facd2 = d2 * vbld_inv(j+nres) + facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres) + DO k = 1, 3 + condor = (erhead_tail(k,2) & + + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j+nres))) gvdwx(k,i) = gvdwx(k,i) & - - dPOLdR1 * hawk - gvdwx(k,j) = gvdwx(k,j) & - + dPOLdR1 * (erhead_tail(k,1) & - -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) + - dPOLdR2 * (erhead_tail(k,2) & + -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) + gvdwx(k,j) = gvdwx(k,j) & + + dPOLdR2 * condor - gvdwc(k,i) = gvdwc(k,i) - dPOLdR1 * erhead_tail(k,1) - gvdwc(k,j) = gvdwc(k,j) + dPOLdR1 * erhead_tail(k,1) + gvdwc(k,i) = gvdwc(k,i) & + - dPOLdR2 * erhead_tail(k,2) + gvdwc(k,j) = gvdwc(k,j) & + + dPOLdR2 * erhead_tail(k,2) END DO - RETURN - END SUBROUTINE eqn - SUBROUTINE enq(Epol) + RETURN + END SUBROUTINE enq + + SUBROUTINE enq_cat(Epol) use calc_data use comm_momo double precision facd3, adler,epol - alphapol2 = alphapol(itypj,itypi) + alphapol2 = alphapolcat(itypj,itypi) !c! R2 - distance between head of jth side chain and tail of ith sidechain R2 = 0.0d0 DO k = 1, 3 @@ -25416,19 +26888,20 @@ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1 !c! dPOLdOM1 = 0.0d0 dPOLdOM2 = 0.0d0 + !c!------------------------------------------------------------------- !c! Return the results !c! (See comments in Eqq) DO k = 1, 3 erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2) END DO - eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) ) + eagle = scalar( erhead_tail(1,2), dC_norm(1,j) ) adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) ) facd2 = d2 * vbld_inv(j+nres) - facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres) + facd3 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres) DO k = 1, 3 condor = (erhead_tail(k,2) & - + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j+nres))) + + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j))) gvdwx(k,i) = gvdwx(k,i) & - dPOLdR2 * (erhead_tail(k,2) & @@ -25443,7 +26916,8 @@ END DO RETURN - END SUBROUTINE enq + END SUBROUTINE enq_cat + SUBROUTINE eqd(Ecl,Elj,Epol) use calc_data use comm_momo @@ -25672,6 +27146,126 @@ END DO RETURN END SUBROUTINE edq + + SUBROUTINE edq_cat(Ecl,Elj,Epol) + use comm_momo + use calc_data + + double precision facd3, adler,ecl,elj,epol + alphapol2 = alphapolcat(itypj,itypi) + w1 = wqdipcat(1,itypi,itypj) + w2 = wqdipcat(2,itypi,itypj) + pis = sig0headcat(itypi,itypj) + eps_head = epsheadcat(itypi,itypj) +!c!------------------------------------------------------------------- +!c! R2 - distance between head of jth side chain and tail of ith sidechain + R2 = 0.0d0 + DO k = 1, 3 +!c! Calculate head-to-tail distances + R2=R2+(chead(k,2)-ctail(k,1))**2 + END DO +!c! Pitagoras + R2 = dsqrt(R2) + +!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj) +!c! & +dhead(1,1,itypi,itypj))**2)) +!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj) +!c! & +dhead(2,1,itypi,itypj))**2)) + + +!c!------------------------------------------------------------------- +!c! ecl + sparrow = w1 * Qi * om1 + hawk = w2 * Qi * Qi * (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 * Qi) / (Rhead**2.0d0) +!c! dF/dom2 + dGCLdOM2 = (2.0d0 * w2 * Qi * Qi * om2) / (Rhead ** 4.0d0) +!c-------------------------------------------------------------------- +!c-------------------------------------------------------------------- +!c Polarization energy +!c Epol + MomoFac2 = (1.0d0 - chi2 * sqom1) + RR2 = R2 * R2 / MomoFac2 + ee2 = exp(-(RR2 / (4.0d0 * a12sq))) + fgb2 = sqrt(RR2 + a12sq * ee2) + epol = 332.0d0 * eps_inout_fac * ((alphapol2/fgb2) ** 4.0d0 ) + dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0) & + / (fgb2 ** 5.0d0) + dFGBdR2 = ( (R2 / MomoFac2) & + * ( 2.0d0 - (0.5d0 * ee2) ) ) & + / (2.0d0 * fgb2) + dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) & + * (2.0d0 - 0.5d0 * ee2) ) & + / (2.0d0 * fgb2) + dPOLdR2 = dPOLdFGB2 * dFGBdR2 +!c! dPOLdR2 = 0.0d0 + dPOLdOM1 = dPOLdFGB2 * dFGBdOM1 +!c! dPOLdOM1 = 0.0d0 + dPOLdOM2 = 0.0d0 +!c!------------------------------------------------------------------- +!c! Elj + pom = (pis / Rhead)**6.0d0 + Elj = 4.0d0 * eps_head * pom * (pom-1.0d0) +!c! derivative of Elj is Glj + dGLJdR = 4.0d0 * eps_head & + * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) & + + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0))) +!c!------------------------------------------------------------------- + +!c! Return the results +!c! (see comments in Eqq) + DO k = 1, 3 + erhead(k) = Rhead_distance(k)/Rhead + erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2) + END DO + erdxi = scalar( erhead(1), dC_norm(1,i+nres) ) + erdxj = scalar( erhead(1), dC_norm(1,j) ) + eagle = scalar( erhead_tail(1,2), dC_norm(1,j) ) + adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) ) + facd1 = d1 * vbld_inv(i+nres) + facd2 = d2 * vbld_inv(j) + facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres) + DO k = 1, 3 + condor = (erhead_tail(k,2) & + + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j))) + + pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres)) + gvdwx(k,i) = gvdwx(k,i) & + - dGCLdR * pom & + - dPOLdR2 * (erhead_tail(k,2) & + -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) & + - dGLJdR * pom + + pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j)) + gvdwx(k,j) = gvdwx(k,j) & + + dGCLdR * pom & + + dPOLdR2 * condor & + + dGLJdR * pom + + + gvdwc(k,i) = gvdwc(k,i) & + - dGCLdR * erhead(k) & + - dPOLdR2 * erhead_tail(k,2) & + - dGLJdR * erhead(k) + + gvdwc(k,j) = gvdwc(k,j) & + + dGCLdR * erhead(k) & + + dPOLdR2 * erhead_tail(k,2) & + + dGLJdR * erhead(k) + + END DO + RETURN + END SUBROUTINE edq_cat + + SUBROUTINE edd(ECL) ! IMPLICIT NONE use comm_momo @@ -25870,4 +27464,178 @@ dPOLdOM2 = 0.0d0 RETURN END SUBROUTINE elgrad_init + + + SUBROUTINE elgrad_init_cat(eheadtail,Egb,Ecl,Elj,Equad,Epol) + use comm_momo + use calc_data + real(kind=8) :: eheadtail,Egb,Ecl,Elj,Equad,Epol,Rb + eps_out=80.0d0 + itypi = itype(i,1) + itypj = itype(j,1) +!c! 1/(Gas Constant * Thermostate temperature) = BetaT +!c! ENABLE THIS LINE WHEN USING CHECKGRAD!!! +!c! t_bath = 300 +!c! BetaT = 1.0d0 / (t_bath * Rb)i + Rb=0.001986d0 + BetaT = 1.0d0 / (298.0d0 * Rb) +!c! Gay-berne var's + sig0ij = sigmacat( itypi,itypj ) + chi1 = chicat( itypi, itypj ) +! chi2 = chi( itypj, itypi ) + chi2 = 0.0d0 +! chi12 = chi1 * chi2 + chi12 = 0.0d0 + chip1 = chippcat( itypi, itypj ) +! chip2 = chipp( itypj, itypi ) + chip2 = 0.0d0 +! chip12 = chip1 * chip2 + chip12 = 0.0d0 +! chi1=0.0 +! chi2=0.0 +! chi12=0.0 +! chip1=0.0 +! chip2=0.0 +! chip12=0.0 +!c! not used by momo potential, but needed by sc_angular which is shared +!c! by all energy_potential subroutines + alf1 = 0.0d0 + alf2 = 0.0d0 + alf12 = 0.0d0 +!c! location, location, location +! xj = c( 1, nres+j ) - xi +! yj = c( 2, nres+j ) - yi +! zj = c( 3, nres+j ) - zi + dxj = dc_norm( 1, nres+j ) + dyj = dc_norm( 2, nres+j ) + dzj = dc_norm( 3, nres+j ) +!c! distance from center of chain(?) to polar/charged head + d1 = dheadcat(1, 1, itypi, itypj) + d2 = dheadcat(2, 1, itypi, itypj) +!c! ai*aj from Fgb + a12sq = rborncat(itypi,itypj) * rborncat(itypj,itypi) +!c! a12sq = a12sq * a12sq +!c! charge of amino acid itypi is... + Qi = ichargecat(itypi) + Qj = ichargecat(itypj) + Qij = Qi * Qj +!c! chis1,2,12 + chis1 = chiscat(itypi,itypj) +! chis2 = chis(itypj,itypi) + chis2 = 0.0d0 +! chis12 = chis1 * chis2 + chis12 = 0.0d0 + sig1 = sigmap1cat(itypi,itypj) + sig2 = sigmap2cat(itypi,itypj) +!c! alpha factors from Fcav/Gcav + b1cav = alphasurcat(1,itypi,itypj) +! b1cav=0.0 + b2cav = alphasurcat(2,itypi,itypj) + b3cav = alphasurcat(3,itypi,itypj) + b4cav = alphasurcat(4,itypi,itypj) + wqd = wquadcat(itypi, itypj) +!c! used by Fgb + eps_in = epsintabcat(itypi,itypj) + eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out)) +!c!------------------------------------------------------------------- +!c! tail location and distance calculations + Rtail = 0.0d0 + DO k = 1, 3 + ctail(k,1)=c(k,i+nres)-dtailcat(1,itypi,itypj)*dc_norm(k,nres+i) + ctail(k,2)=c(k,j+nres)-dtailcat(2,itypi,itypj)*dc_norm(k,nres+j) + END DO +!c! tail distances will be themselves usefull elswhere +!c1 (in Gcav, for example) + Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 ) + Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 ) + Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 ) + Rtail = dsqrt( & + (Rtail_distance(1)*Rtail_distance(1)) & + + (Rtail_distance(2)*Rtail_distance(2)) & + + (Rtail_distance(3)*Rtail_distance(3))) +!c!------------------------------------------------------------------- +!c! Calculate location and distance between polar heads +!c! distance between heads +!c! for each one of our three dimensional space... + d1 = dheadcat(1, 1, itypi, itypj) + d2 = dheadcat(2, 1, itypi, itypj) + + DO k = 1,3 +!c! location of polar head is computed by taking hydrophobic centre +!c! and moving by a d1 * dc_norm vector +!c! see unres publications for very informative images + chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres) + chead(k,2) = c(k, j+nres) + d2 * dc_norm(k, j+nres) +!c! distance +!c! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres)) +!c! Rsc(k) = Rsc_distance(k) * Rsc_distance(k) + Rhead_distance(k) = chead(k,2) - chead(k,1) + END DO +!c! pitagoras (root of sum of squares) + Rhead = dsqrt( & + (Rhead_distance(1)*Rhead_distance(1)) & + + (Rhead_distance(2)*Rhead_distance(2)) & + + (Rhead_distance(3)*Rhead_distance(3))) +!c!------------------------------------------------------------------- +!c! zero everything that should be zero'ed + Egb = 0.0d0 + ECL = 0.0d0 + Elj = 0.0d0 + Equad = 0.0d0 + Epol = 0.0d0 + eheadtail = 0.0d0 + dGCLdOM1 = 0.0d0 + dGCLdOM2 = 0.0d0 + dGCLdOM12 = 0.0d0 + dPOLdOM1 = 0.0d0 + dPOLdOM2 = 0.0d0 + RETURN + END SUBROUTINE elgrad_init_cat + + + double precision function tschebyshev(m,n,x,y) + implicit none + integer i,m,n + double precision x(n),y,yy(0:maxvar),aux +!c Tschebyshev polynomial. Note that the first term is omitted +!c m=0: the constant term is included +!c m=1: the constant term is not included + yy(0)=1.0d0 + yy(1)=y + do i=2,n + yy(i)=2*yy(1)*yy(i-1)-yy(i-2) + enddo + aux=0.0d0 + do i=m,n + aux=aux+x(i)*yy(i) + enddo + tschebyshev=aux + return + end function tschebyshev +!C-------------------------------------------------------------------------- + double precision function gradtschebyshev(m,n,x,y) + implicit none + integer i,m,n + double precision x(n+1),y,yy(0:maxvar),aux +!c Tschebyshev polynomial. Note that the first term is omitted +!c m=0: the constant term is included +!c m=1: the constant term is not included + yy(0)=1.0d0 + yy(1)=2.0d0*y + do i=2,n + yy(i)=2*y*yy(i-1)-yy(i-2) + enddo + aux=0.0d0 + do i=m,n + aux=aux+x(i+1)*yy(i)*(i+1) +!C print *, x(i+1),yy(i),i + enddo + gradtschebyshev=aux + return + end function gradtschebyshev + + + + + end module energy