- module energy
+ module energy
!-----------------------------------------------------------------------------
use io_units
use names
! 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
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
! include 'COMMON.TIME1'
real(kind=8) :: time00
!el local variables
- integer :: n_corr,n_corr1,ierror
+ integer :: n_corr,n_corr1,ierror,imatupdate
real(kind=8) :: etors,edihcnstr,etors_d,esccor,ehpb
real(kind=8) :: evdw,evdw1,evdw2,evdw2_14,escloc,ees,eel_loc
real(kind=8) :: eello_turn3,eello_turn4,estr,ebe,eliptran,etube, &
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
integer ishield_listbuf(-1:nres), &
shield_listbuf(maxcontsshi,-1:nres),k,j,i,iii,impishi,mojint,jjj
-
-
+! print *,"I START ENERGY"
+ imatupdate=100
+! if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
! real(kind=8), dimension(:),allocatable:: fac_shieldbuf
! real(kind=8), dimension(:,:,:),allocatable:: &
! grad_shield_locbuf,grad_shield_sidebuf
! 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
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)
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
! call chainbuild_cart
endif
+! print *,"itime_mat",itime_mat,imatupdate
+ if (nfgtasks.gt.1) then
+ call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
+ endif
+ if (mod(itime_mat,imatupdate).eq.0) call make_SCp_inter_list
+ if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
+ if (mod(itime_mat,imatupdate).eq.0) call make_pp_inter_list
+
! print *,'Processor',myrank,' calling etotal ipot=',ipot
! print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
#else
.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
! 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"
! 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"
!
call AFMforce(Eafmforce)
else if (selfguide.gt.0) then
call AFMvel(Eafmforce)
+ else
+ Eafmforce=0.0d0
endif
endif
if (tubemode.eq.1) then
etors_nucl=0.0d0
estr_nucl=0.0d0
ecorr3_nucl=0.0d0
+ ecorr_nucl=0.0d0
ebe_nucl=0.0d0
evdwsb=0.0d0
eelsb=0.0d0
eelpsb=0.0d0
evdwpp=0.0d0
eespp=0.0d0
+ etors_d_nucl=0.0d0
endif
! write(iout,*) ecorr_nucl,"ecorr_nucl",nres_molec(2)
+! print *,"before ecatcat",wcatcat
+ if (nres_molec(5).gt.0) then
if (nfgtasks.gt.1) then
if (fg_rank.eq.0) then
call ecatcat(ecationcation)
else
call ecatcat(ecationcation)
endif
+ if (oldion.gt.0) then
call ecat_prot(ecation_prot)
- if (nres_molec(2).gt.0) then
+ else
+ call ecats_prot_amber(ecation_prot)
+ endif
+ else
+ ecationcation=0.0d0
+ ecation_prot=0.0d0
+ endif
+ if ((nres_molec(2).gt.0).and.(nres_molec(1).gt.0)) then
call eprot_sc_base(escbase)
call epep_sc_base(epepbase)
call eprot_sc_phosphate(escpho)
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
! 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"
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
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
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
!-----------------------------------------------------------------------------
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)
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,&
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,&
! include 'COMMON.SBRIDGE'
logical :: lprn
!el local variables
- integer :: iint,itypi,itypi1,itypj,subchap
+ integer :: iint,itypi,itypi1,itypj,subchap,icont
real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
real(kind=8) :: evdw,sig0ij
real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
dPOLdOM1=0.0d0
- do i=iatsc_s,iatsc_e
+ do icont=g_listscsc_start,g_listscsc_end
+ i=newcontlisti(icont)
+ j=newcontlistj(icont)
+
+! do i=iatsc_s,iatsc_e
!C print *,"I am in EVDW",i
itypi=iabs(itype(i,1))
! if (i.ne.47) cycle
!
! Calculate SC interaction energy.
!
- do iint=1,nint_gr(i)
- do j=istart(i,iint),iend(i,iint)
+! do iint=1,nint_gr(i)
+! do j=istart(i,iint),iend(i,iint)
IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
call dyn_ssbond_ene(i,j,evdwij)
evdw=evdw+evdwij
! 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
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
! Calculate angular part of the gradient.
call sc_grad
ENDIF ! dyn_ss
- enddo ! j
- enddo ! iint
+! enddo ! j
+! enddo ! iint
enddo ! i
! print *,"ZALAMKA", evdw
! write (iout,*) "Number of loop steps in EGB:",ind
! 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),nloctyp
+! print *, "i,",molnum(i),i,itype(i-2,1)
+ if (molnum(i).eq.1) then
+ if (itype(i-2,1).eq.ntyp1) then
+ iti=nloctyp
+ else
+ iti = itype2loc(itype(i-2,1))
+ endif
+ 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))
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
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))
0.0d0,1.0d0,0.0d0,&
0.0d0,0.0d0,1.0d0/),shape(unmat))
!el local variables
- integer :: i,k,j
+ integer :: i,k,j,icont
real(kind=8) :: ees,evdw1,eel_loc,eello_turn3,eello_turn4
real(kind=8) :: fac,t_eelecij,fracinbuf
! Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
!
! print *,"iatel_s,iatel_e,",iatel_s,iatel_e
- do i=iatel_s,iatel_e
+! do i=iatel_s,iatel_e
+! JPRDLC
+ do icont=g_listpp_start,g_listpp_end
+ i=newcontlistppi(icont)
+ j=newcontlistppj(icont)
if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
- do j=ielstart(i),ielend(i)
+! do j=ielstart(i),ielend(i)
! write (iout,*) i,j,itype(i,1),itype(j,1)
if (itype(j,1).eq.ntyp1.or. itype(j+1,1).eq.ntyp1) cycle
call eelecij(i,j,ees,evdw1,eel_loc)
- enddo ! j
+! enddo ! j
num_cont_hb(i)=num_conti
enddo ! i
! write (iout,*) "Number of loop steps in EELEC:",ind
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
! sss_ele_grad=0.0d0
! print *,sss_ele_cut,sss_ele_grad,&
! (rij),r_cut_ele,rlamb_ele
-! if (sss_ele_cut.le.0.0) go to 128
+ if (sss_ele_cut.le.0.0) go to 128
rmij=1.0D0/rij
r3ij=rrmij*rmij
!grad enddo
!grad enddo
! 9/28/08 AL Gradient compotents will be summed only at the end
- ggg(1)=facvdw*xj &
+ ggg(1)=facvdw*xj+sss_ele_grad*rmij*evdwij*xj &
*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
- ggg(2)=facvdw*yj &
+ ggg(2)=facvdw*yj+sss_ele_grad*rmij*evdwij*yj &
*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
- ggg(3)=facvdw*zj &
+ ggg(3)=facvdw*zj+sss_ele_grad*rmij*evdwij*zj &
*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
do k=1,3
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
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 &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+
+!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 &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+
+!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 &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+
+ 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 &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+#endif
! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
! eel_loc_ij=0.0
ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij) &
*sss_ele_cut &
*fac_shield(i)*fac_shield(j)
+! *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) &
*sss_ele_cut &
*fac_shield(i)*fac_shield(j)
+! *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
! Diagnostics. Comment out or remove after debugging!
! ees0p(num_conti,i)=0.5D0*fac3*ees0pij
gacontp_hb1(k,num_conti,i)= & !ghalfp+
(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i)) &
+ ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) &
- *sss_ele_cut*fac_shield(i)*fac_shield(j)
+ *sss_ele_cut*fac_shield(i)*fac_shield(j) ! &
+! *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
gacontp_hb2(k,num_conti,i)= & !ghalfp+
(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j)) &
+ ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)&
- *sss_ele_cut*fac_shield(i)*fac_shield(j)
+ *sss_ele_cut*fac_shield(i)*fac_shield(j)! &
+! *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
gacontp_hb3(k,num_conti,i)=gggp(k) &
*sss_ele_cut*fac_shield(i)*fac_shield(j)
+! *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
gacontm_hb1(k,num_conti,i)= & !ghalfm+
(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i)) &
+ ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) &
*sss_ele_cut*fac_shield(i)*fac_shield(j)
+! *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
gacontm_hb2(k,num_conti,i)= & !ghalfm+
(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j)) &
+ ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) &
*sss_ele_cut*fac_shield(i)*fac_shield(j)
+! *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
gacontm_hb3(k,num_conti,i)=gggm(k) &
*sss_ele_cut*fac_shield(i)*fac_shield(j)
+! *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
enddo
! Diagnostics. Comment out or remove after debugging!
! 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
!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
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) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+ 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) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+
+!C#endif
+
+
+
if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. &
(shield_mode.gt.0)) then
!C print *,i,j
! 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,&
!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
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
! 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),
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) &
! include 'COMMON.CONTROL'
real(kind=8),dimension(3) :: ggg
!el local variables
- integer :: i,iint,j,k,iteli,itypj,subchap
+ integer :: i,iint,j,k,iteli,itypj,subchap,icont
real(kind=8) :: evdw2,evdw2_14,xi,yi,zi,xj,yj,zj,rrij,fac,&
e1,e2,evdwij,rij
real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
evdw2_14=0.0d0
!d print '(a)','Enter ESCP'
!d write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
- do i=iatscp_s,iatscp_e
+! do i=iatscp_s,iatscp_e
+ do icont=g_listscp_start,g_listscp_end
+ i=newcontlistscpi(icont)
+ j=newcontlistscpj(icont)
if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
zi=mod(zi,boxzsize)
if (zi.lt.0) zi=zi+boxzsize
- do iint=1,nscp_gr(i)
+! do iint=1,nscp_gr(i)
- do j=iscpstart(i,iint),iscpend(i,iint)
+! do j=iscpstart(i,iint),iscpend(i,iint)
itypj=iabs(itype(j,1))
if (itypj.eq.ntyp1) cycle
! Uncomment following three lines for SC-p interactions
gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
enddo
- enddo
+! enddo
- enddo ! iint
+! enddo ! iint
enddo ! i
do i=1,nct
do j=1,3
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.
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
! & dscp1,dscp2,sumene
! sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
escloc = escloc + sumene
+ if (energy_dec) write (2,*) "i",i," itype",itype(i,1)," it",it, &
+ " escloc",sumene,escloc,it,itype(i,1)
! write (2,*) "i",i," escloc",sumene,escloc,it,itype(i,1)
! & ,zz,xx,yy
!#define DEBUG
end subroutine etor_d
#else
!-----------------------------------------------------------------------------
- subroutine etor(etors,edihcnstr)
+ subroutine etor(etors)
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! include 'COMMON.VAR'
! 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)
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
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
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))
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)
#ifdef TIMING
time01=MPI_Wtime()
#endif
+!#define DEBUG
#ifdef DEBUG
write (iout,*) "sum_gradient gvdwc, gvdwx"
do i=1,nres
! 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
+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)
+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) &
enddo
#endif
!#undef DEBUG
- do i=1,nres
+ do i=0,nres
do j=1,3
gloc_scbuf(j,i)=gloc_sc(j,i,icg)
enddo
call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres,&
MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
time_reduce=time_reduce+MPI_Wtime()-time00
- call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres,&
+ call MPI_Reduce(gloc_scbuf(1,0),gloc_sc(1,0,icg),3*nres+3,&
MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
time_reduce=time_reduce+MPI_Wtime()-time00
!#define DEBUG
! print *,"gradbuf",gradbufc(1,1),gradc(1,1,icg)
#ifdef DEBUG
write (iout,*) "gloc_sc after reduce"
- do i=1,nres
+ do i=0,nres
do j=1,1
write (iout,*) i,j,gloc_sc(j,i,icg)
enddo
enddo
return
end subroutine sc_grad
-#ifdef CRYST_THETA
-!-----------------------------------------------------------------------------
- subroutine mixder(thetai,thet_pred_mean,theta0i,E_tc_t)
- use comm_calcthet
-! implicit real*8 (a-h,o-z)
+ subroutine sc_grad_cat
+ use calc_data
+ real(kind=8), dimension(3) :: dcosom1,dcosom2
+ eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 &
+ +dCAVdOM1+ dGCLdOM1+ dPOLdOM1
+ eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 &
+ +dCAVdOM2+ dGCLdOM2+ dPOLdOM2
+
+ eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 &
+ -2.0D0*alf12*eps3der+sigder*sigsq_om12&
+ +dCAVdOM12+ dGCLdOM12
+! diagnostics only
+! eom1=0.0d0
+! eom2=0.0d0
+! eom12=evdwij*eps1_om12
+! end diagnostics
+
+ do k=1,3
+ dcosom1(k)=rij*(dc_norm(k,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
+ gradpepcatx(k,i)=gradpepcatx(k,i)-gg(k) &
+ +(eom12*(dc_norm(k,j)-om12*dc_norm(k,nres+i)) &
+ +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+
+! gradpepcatx(k,j)=gradpepcatx(k,j)+gg(k) &
+! +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,j)) &
+! +eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv
+
+! write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
+! +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+! 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
+!
+ do l=1,3
+ gradpepcat(l,i)=gradpepcat(l,i)-gg(l)
+ gradpepcat(l,j)=gradpepcat(l,j)+gg(l)
+ enddo
+ end subroutine sc_grad_cat
+
+ subroutine sc_grad_cat_pep
+ use calc_data
+ real(kind=8), dimension(3) :: dcosom1,dcosom2
+ eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 &
+ +dCAVdOM1+ dGCLdOM1+ dPOLdOM1
+ eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 &
+ +dCAVdOM2+ dGCLdOM2+ dPOLdOM2
+
+ eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 &
+ -2.0D0*alf12*eps3der+sigder*sigsq_om12&
+ +dCAVdOM12+ dGCLdOM12
+! diagnostics only
+! eom1=0.0d0
+! eom2=0.0d0
+! eom12=evdwij*eps1_om12
+! end diagnostics
+
+ do k=1,3
+ dcosom1(k) = rij * (dc_norm(k,i) - om1 * erij(k))
+ dcosom2(k) = rij * (dc_norm(k,nres+j) - om2 * erij(k))
+ gg(k) = gg(k) + eom1 * dcosom1(k) + eom2 * dcosom2(k)
+ gvdwc_pepbase(k,i)= gvdwc_pepbase(k,i) +0.5*(- gg(k)) &
+ + (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i)))&
+ *dsci_inv*2.0 &
+ - (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0
+ gvdwc_pepbase(k,i+1)= gvdwc_pepbase(k,i+1) +0.5*(- gg(k)) &
+ - (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i))) &
+ *dsci_inv*2.0 &
+ + (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0
+ gradpepcat(k,j)=gradpepcat(k,j)+gg(k)
+ enddo
+ end subroutine sc_grad_cat_pep
+
+#ifdef CRYST_THETA
+!-----------------------------------------------------------------------------
+ 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'
! call intcartderiv
! call checkintcartgrad
call zerograd
- aincr=1.0D-4
+ aincr=1.0D-5
write(iout,*) 'Calling CHECK_ECARTINT.'
nf=0
icall=0
! call intcartderiv
! call checkintcartgrad
call zerograd
- aincr=2.0D-5
+ aincr=1.0D-6
write(iout,*) 'Calling CHECK_ECARTINT.',aincr
nf=0
icall=0
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
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
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
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
!
! 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
!
!#define DEBUG
!el write (iout,*) "After sum_gradient"
#ifdef DEBUG
-!el write (iout,*) "After sum_gradient"
+ write (iout,*) "After sum_gradient"
do i=1,nres-1
write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3)
write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3)
dphi(j,1,i)=0.0d0
dphi(j,2,i)=0.0d0
dphi(j,3,i)=0.0d0
+ dcosomicron(j,1,1,i)=0.0d0
+ dcosomicron(j,1,2,i)=0.0d0
+ dcosomicron(j,2,1,i)=0.0d0
+ dcosomicron(j,2,2,i)=0.0d0
enddo
enddo
! Derivatives of theta's
#else
do i=3,nres
#endif
- if ((itype(i-1,1).ne.10).and.(itype(i-1,1).ne.ntyp1)) then
+ if ((itype(i-1,1).ne.10).and.(itype(i-1,1).ne.ntyp1).and.molnum(i).ne.5) then
cost1=dcos(omicron(1,i))
sint1=sqrt(1-cost1*cost1)
cost2=dcos(omicron(2,i))
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
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))
allocate(uygrad(3,3,2,nres))
allocate(uzgrad(3,3,2,nres))
!(3,3,2,maxres)
+! allocateion of lists JPRDLA
+ allocate(newcontlistppi(200*nres))
+ allocate(newcontlistscpi(200*nres))
+ allocate(newcontlisti(200*nres))
+ allocate(newcontlistppj(200*nres))
+ allocate(newcontlistscpj(200*nres))
+ allocate(newcontlistj(200*nres))
return
end subroutine alloc_ener_arrays
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
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,
!c------------------------------------------------------------------------------
#endif
subroutine ecatcat(ecationcation)
- integer :: i,j,itmp,xshift,yshift,zshift,subchap,k
+ integer :: i,j,itmp,xshift,yshift,zshift,subchap,k,itypi,itypj
real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,&
r7,r4,ecationcation,k0,rcal
real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
epscalc=0.05
r06 = rcat0**6
r012 = r06**2
- k0 = 332.0*(2.0*2.0)/80.0
+! k0 = 332.0*(2.0*2.0)/80.0
itmp=0
do i=1,4
xi=c(1,i)
yi=c(2,i)
zi=c(3,i)
-
+! write (iout,*) i,"TUTUT",c(1,i)
+ itypi=itype(i,5)
xi=mod(xi,boxxsize)
if (xi.lt.0) xi=xi+boxxsize
yi=mod(yi,boxysize)
if (zi.lt.0) zi=zi+boxzsize
do j=i+1,itmp+nres_molec(5)
+ itypj=itype(j,5)
+! print *,i,j,itypi,itypj
+ k0 = 332.0*(ichargecat(itypi)*ichargecat(itypj))/80.0
! print *,i,j,'catcat'
xj=c(1,j)
yj=c(2,j)
! r06 = rcat0**6
! r012 = r06**2
! k0 = 332*(2*2)/80
- Evan1cat=epscalc*(r012/rcal**6)
- Evan2cat=epscalc*2*(r06/rcal**3)
+ Evan1cat=epscalc*(r012/(rcal**6))
+ Evan2cat=epscalc*2*(r06/(rcal**3))
Eeleccat=k0/ract
r7 = rcal**7
r4 = rcal**4
gradcatcat(k,i)=gradcatcat(k,i)-gg(k)
gradcatcat(k,j)=gradcatcat(k,j)+gg(k)
enddo
-
+ if (energy_dec) write (iout,*) i,j,Evan1cat,Evan2cat,Eeleccat,&
+ r012,rcal**6,ichargecat(itypi)*ichargecat(itypj)
! write(iout,*) "ecatcat",i,j, ecationcation,xj,yj,zj
ecationcation=ecationcation+Evan1cat+Evan2cat+Eeleccat
enddo
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(evdw)
+! 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)
+
+ evdw=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
+! go to 17
! 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).or.(itypi.eq.10)) 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,5))
+ if ((itypj.eq.ntyp1)) cycle
+ CALL elgrad_init_cat(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+
+ dscj_inv=0.0
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
- xj=dmod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=dmod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=dmod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
xj=xj_safe+xshift*boxxsize
yj=yj_safe+yshift*boxysize
zj=zj_safe+zshift*boxzsize
zj_temp=zj
subchap=1
endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
xj=xj_temp-xi
yj=yj_temp-yi
zj=zj_temp-zi
- else
+ else
xj=xj_safe-xi
yj=yj_safe-yi
zj=zj_safe-zi
- endif
-! enddo
-! enddo
- rcpm = sqrt(xj**2+yj**2+zj**2)
- drcp_norm(1)=xj/rcpm
- drcp_norm(2)=yj/rcpm
- drcp_norm(3)=zj/rcpm
- dcmag=0.0
- do k=1,3
- dcmag=dcmag+dc(k,i)**2
- enddo
- dcmag=dsqrt(dcmag)
- do k=1,3
- myd_norm(k)=dc(k,i)/dcmag
- enddo
- costhet=drcp_norm(1)*myd_norm(1)+drcp_norm(2)*myd_norm(2)+&
- drcp_norm(3)*myd_norm(3)
- rsecp = rcpm**2
- Ir = 1.0d0/rcpm
- Irsecp = 1.0d0/rsecp
- Irthrp = Irsecp/rcpm
- Irfourp = Irthrp/rcpm
- Irfiftp = Irfourp/rcpm
- Irsistp=Irfiftp/rcpm
- Irseven=Irsistp/rcpm
- Irtwelv=Irsistp*Irsistp
- Irthir=Irtwelv/rcpm
+ 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 = chi1cat(itypi,itypj)
+ chis1 = chis1cat(itypi,itypj)
+ chip1 = chipp1cat(itypi,itypj)
+! chi1=0.0d0
+! chis1=0.0d0
+! chip1=0.0d0
+ chi2=0.0
+ chip2=0.0
+ chis2=0.0
+! chis2 = chis(itypj,itypi)
+ chis12 = chis1 * chis2
+ sig1 = sigmap1cat(itypi,itypj)
+! sig2 = sigmap2(itypi,itypj)
+! alpha factors from Fcav/Gcav
+ b1cav = alphasurcat(1,itypi,itypj)
+ b2cav = alphasurcat(2,itypi,itypj)
+ b3cav = alphasurcat(3,itypi,itypj)
+ b4cav = alphasurcat(4,itypi,itypj)
+
+! used to determine whether we want to do quadrupole calculations
+ eps_in = epsintabcat(itypi,itypj)
+ if (eps_in.eq.0.0) eps_in=1.0
+
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+! Rtail = 0.0d0
+
+ DO k = 1, 3
+ ctail(k,1)=c(k,i+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_cat(itypi,itypj)
+! print *,"ADAM",aa_aq(itypi,itypj)
+
+! c1 = 0.0d0
+ c2 = fac * bb_aq_cat(itypi,itypj)
+! c2 = 0.0d0
+ evdwij = eps1 * eps2rt * eps3rt * ( c1 + c2 )
+ eps2der = eps3rt * evdwij
+ eps3der = eps2rt * evdwij
+! evdwij = 4.0d0 * eps2rt * eps3rt * evdwij
+ evdwij = eps2rt * eps3rt * evdwij
+!#ifdef TSCSC
+! IF (bb_aq(itypi,itypj).gt.0) THEN
+! evdw_p = evdw_p + evdwij
+! ELSE
+! evdw_m = evdw_m + evdwij
+! END IF
+!#else
+ evdw = evdw &
+ + evdwij
+!#endif
+ c1 = c1 * eps1 * eps2rt**2 * eps3rt**2
+ fac = -expon * (c1 + evdwij) * rij_shift
+ sigder = fac * sigder
+! Calculate distance derivative
+ gg(1) = fac
+ gg(2) = fac
+ gg(3) = fac
+
+ fac = chis1 * sqom1 + chis2 * sqom2 &
+ - 2.0d0 * chis12 * om1 * om2 * om12
+ pom = 1.0d0 - chis1 * chis2 * sqom12
+ 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 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres)
+ facd2 = dtailcat(2,itypi,itypj) * vbld_inv(j+nres)
+ DO k = 1, 3
+ pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres))
+ gradpepcatx(k,i) = gradpepcatx(k,i) &
+ - (( dFdR + gg(k) ) * pom)
+ pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres))
+! gvdwx(k,j) = gvdwx(k,j) &
+! + (( dFdR + gg(k) ) * pom)
+ gradpepcat(k,i) = gradpepcat(k,i) &
+ - (( dFdR + gg(k) ) * ertail(k))
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ + (( dFdR + gg(k) ) * ertail(k))
+ gg(k) = 0.0d0
+ ENDDO
+!c! Compute head-head and head-tail energies for each state
+ isel = iabs(Qi) + 1 ! ion is always charged so iabs(Qj)
+ IF (isel.eq.0) THEN
+!c! No charges - do nothing
+ eheadtail = 0.0d0
+
+ ELSE IF (isel.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) 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
+ write(iout,*) "KURWA0",d1
+
+ CALL edq_cat(ecl, elj, epol)
+ eheadtail = ECL + elj + epol
+! eheadtail = 0.0d0
+
+ ELSE IF ((isel.eq.2)) 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 ! i
+!c write (iout,*) "Number of loop steps in EGB:",ind
+!c energy_dec=.false.
+! print *,"EVDW KURW",evdw,nres
+!!! return
+ 17 continue
+ do i=ibond_start,ibond_end
+
+! print *,"I am in EVDW",i
+ itypi=10 ! the peptide group parameters are for glicine
+
+! if (i.ne.47) cycle
+ if ((itype(i,1).eq.ntyp1).or.itype(i+1,1).eq.ntyp1) cycle
+ itypi1=iabs(itype(i+1,1))
+ xi=(c(1,i)+c(1,i+1))/2.0
+ yi=(c(2,i)+c(2,i+1))/2.0
+ zi=(c(3,i)+c(3,i+1))/2.0
+ xi=dmod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=dmod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=dmod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+ dxi=dc_norm(1,i)
+ dyi=dc_norm(2,i)
+ dzi=dc_norm(3,i)
+ dsci_inv=vbld_inv(i+1)/2.0
+ do j=itmp+1,itmp+nres_molec(5)
+
+! Calculate SC interaction energy.
+ itypj=iabs(itype(j,5))
+ if ((itypj.eq.ntyp1)) cycle
+ CALL elgrad_init_cat_pep(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+
+ dscj_inv=0.0
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+
+ dxj = 0.0d0! dc_norm( 1, nres+j )
+ dyj = 0.0d0!dc_norm( 2, nres+j )
+ dzj = 0.0d0! dc_norm( 3, nres+j )
+
+ itypi = 10
+ itypj = itype(j,5)
+! Parameters from fitting the analitical expressions to the PMF obtained by umbrella
+! sampling performed with amber package
+! alf1 = 0.0d0
+! alf2 = 0.0d0
+! alf12 = 0.0d0
+! a12sq = rborn(itypi,itypj) * rborn(itypj,itypi)
+ chi1 = chi1cat(itypi,itypj)
+ chis1 = chis1cat(itypi,itypj)
+ chip1 = chipp1cat(itypi,itypj)
+! chi1=0.0d0
+! chis1=0.0d0
+! chip1=0.0d0
+ chi2=0.0
+ chip2=0.0
+ chis2=0.0
+! chis2 = chis(itypj,itypi)
+ chis12 = chis1 * chis2
+ sig1 = sigmap1cat(itypi,itypj)
+! sig2 = sigmap2(itypi,itypj)
+! alpha factors from Fcav/Gcav
+ b1cav = alphasurcat(1,itypi,itypj)
+ b2cav = alphasurcat(2,itypi,itypj)
+ b3cav = alphasurcat(3,itypi,itypj)
+ b4cav = alphasurcat(4,itypi,itypj)
+
+! used to determine whether we want to do quadrupole calculations
+ eps_in = epsintabcat(itypi,itypj)
+ if (eps_in.eq.0.0) eps_in=1.0
+
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+! Rtail = 0.0d0
+
+ DO k = 1, 3
+ ctail(k,1)=(c(k,i)+c(k,i+1))/2.0
+ ctail(k,2)=c(k,j)
+ END DO
+!c! tail distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+ Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 )
+ Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 )
+ Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 )
+ Rtail = dsqrt( &
+ (Rtail_distance(1)*Rtail_distance(1)) &
+ + (Rtail_distance(2)*Rtail_distance(2)) &
+ + (Rtail_distance(3)*Rtail_distance(3)))
+! tail location and distance calculations
+! dhead1
+ d1 = dheadcat(1, 1, itypi, itypj)
+! print *,"d1",d1
+! d1=0.0d0
+! d2 = dhead(2, 1, itypi, itypj)
+ DO k = 1,3
+! location of polar head is computed by taking hydrophobic centre
+! and moving by a d1 * dc_norm vector
+! see unres publications for very informative images
+ chead(k,1) = (c(k, i)+c(k,i+1))/2.0 + d1 * dc_norm(k, i)
+ chead(k,2) = c(k, j)
+! distance
+! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ Rhead_distance(k) = chead(k,2) - chead(k,1)
+ END DO
+! pitagoras (root of sum of squares)
+ Rhead = dsqrt( &
+ (Rhead_distance(1)*Rhead_distance(1)) &
+ + (Rhead_distance(2)*Rhead_distance(2)) &
+ + (Rhead_distance(3)*Rhead_distance(3)))
+!-------------------------------------------------------------------
+! zero everything that should be zero'ed
+ evdwij = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ Fcav=0.0d0
+ eheadtail = 0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ Fcav = 0.0d0
+ dFdR = 0.0d0
+ dCAVdOM1 = 0.0d0
+ dCAVdOM2 = 0.0d0
+ dCAVdOM12 = 0.0d0
+ dscj_inv = vbld_inv(j+nres)
+! print *,i,j,dscj_inv,dsci_inv
+! rij holds 1/(distance of Calpha atoms)
+ rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj)
+ rij = dsqrt(rrij)
+ CALL sc_angular
+! this should be in elgrad_init but om's are calculated by sc_angular
+! which in turn is used by older potentials
+! om = omega, sqom = om^2
+ sqom1 = om1 * om1
+ sqom2 = om2 * om2
+ sqom12 = om12 * om12
+
+! now we calculate EGB - Gey-Berne
+! It will be summed up in evdwij and saved in evdw
+ sigsq = 1.0D0 / sigsq
+ sig = sig0ij * dsqrt(sigsq)
+! rij_shift = 1.0D0 / rij - sig + sig0ij
+ rij_shift = Rtail - sig + sig0ij
+ IF (rij_shift.le.0.0D0) THEN
+ evdw = 1.0D20
+ RETURN
+ END IF
+ sigder = -sig * sigsq
+ rij_shift = 1.0D0 / rij_shift
+ fac = rij_shift**expon
+ c1 = fac * fac * aa_aq_cat(itypi,itypj)
+! print *,"ADAM",aa_aq(itypi,itypj)
+
+! c1 = 0.0d0
+ c2 = fac * bb_aq_cat(itypi,itypj)
+! c2 = 0.0d0
+ evdwij = eps1 * eps2rt * eps3rt * ( c1 + c2 )
+ eps2der = eps3rt * evdwij
+ eps3der = eps2rt * evdwij
+! evdwij = 4.0d0 * eps2rt * eps3rt * evdwij
+ evdwij = eps2rt * eps3rt * evdwij
+!#ifdef TSCSC
+! IF (bb_aq(itypi,itypj).gt.0) THEN
+! evdw_p = evdw_p + evdwij
+! ELSE
+! evdw_m = evdw_m + evdwij
+! END IF
+!#else
+ evdw = evdw &
+ + evdwij
+!#endif
+ c1 = c1 * eps1 * eps2rt**2 * eps3rt**2
+ fac = -expon * (c1 + evdwij) * rij_shift
+ sigder = fac * sigder
+! Calculate distance derivative
+ gg(1) = fac
+ gg(2) = fac
+ gg(3) = fac
+
+ fac = chis1 * sqom1 + chis2 * sqom2 &
+ - 2.0d0 * chis12 * om1 * om2 * om12
+
+ pom = 1.0d0 - chis1 * chis2 * sqom12
+! print *,"TUT2",fac,chis1,sqom1,pom
+ Lambf = (1.0d0 - (fac / pom))
+ Lambf = dsqrt(Lambf)
+ sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0)
+ Chif = Rtail * sparrow
+ ChiLambf = Chif * Lambf
+ eagle = dsqrt(ChiLambf)
+ bat = ChiLambf ** 11.0d0
+ top = b1cav * ( eagle + b2cav * ChiLambf - b3cav )
+ bot = 1.0d0 + b4cav * (ChiLambf ** 12.0d0)
+ botsq = bot * bot
+ Fcav = top / bot
+
+ dtop = b1cav * ((Lambf / (2.0d0 * eagle)) + (b2cav * Lambf))
+ dbot = 12.0d0 * b4cav * bat * Lambf
+ dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow
+
+ dtop = b1cav * ((Chif / (2.0d0 * eagle)) + (b2cav * Chif))
+ dbot = 12.0d0 * b4cav * bat * Chif
+ eagle = Lambf * pom
+ dFdOM1 = -(chis1 * om1 - chis12 * om2 * om12) / (eagle)
+ dFdOM2 = -(chis2 * om2 - chis12 * om1 * om12) / (eagle)
+ dFdOM12 = chis12 * (chis1 * om1 * om12 - om2) &
+ * (chis2 * om2 * om12 - om1) / (eagle * pom)
+
+ dFdL = ((dtop * bot - top * dbot) / botsq)
+ dCAVdOM1 = dFdL * ( dFdOM1 )
+ dCAVdOM2 = dFdL * ( dFdOM2 )
+ dCAVdOM12 = dFdL * ( dFdOM12 )
+
+ DO k= 1, 3
+ ertail(k) = Rtail_distance(k)/Rtail
+ END DO
+ erdxi = scalar( ertail(1), dC_norm(1,i) )
+ erdxj = scalar( ertail(1), dC_norm(1,j) )
+ facd1 = dtailcat(1,itypi,itypj) * vbld_inv(i)
+ facd2 = dtailcat(2,itypi,itypj) * vbld_inv(j+nres)
+ DO k = 1, 3
+ pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i))
+! gradpepcatx(k,i) = gradpepcatx(k,i) &
+! - (( dFdR + gg(k) ) * pom)
+ pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres))
+! gvdwx(k,j) = gvdwx(k,j) &
+! + (( dFdR + gg(k) ) * pom)
+ gradpepcat(k,i) = gradpepcat(k,i) &
+ - (( dFdR + gg(k) ) * ertail(k))/2.0d0
+ gradpepcat(k,i+1) = gradpepcat(k,i+1) &
+ - (( dFdR + gg(k) ) * ertail(k))/2.0d0
+
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ + (( dFdR + gg(k) ) * ertail(k))
+ gg(k) = 0.0d0
+ ENDDO
+!c! Compute head-head and head-tail energies for each state
+ isel = 3
+!c! Dipole-charge interactions
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ Qi=Qi*2
+ Qij=Qij*2
+ endif
+ if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+ Qj=Qj*2
+ Qij=Qij*2
+ endif
+ CALL edq_cat_pep(ecl, elj, epol)
+ eheadtail = ECL + elj + epol
+! print *,"i,",i,eheadtail
+! eheadtail = 0.0d0
+
+ evdw = evdw + Fcav + eheadtail
+
+ IF (energy_dec) write (iout,'(2(1x,a3,i3),3f6.2,10f16.7)') &
+ restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
+ 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,&
+ Equad,evdwij+Fcav+eheadtail,evdw
+! evdw = evdw + Fcav + eheadtail
+
+! iF (nstate(itypi,itypj).eq.1) THEN
+ CALL sc_grad_cat_pep
+! END IF
+!c!-------------------------------------------------------------------
+!c! NAPISY KONCOWE
+ END DO ! j
+ END DO ! i
+!c write (iout,*) "Number of loop steps in EGB:",ind
+!c energy_dec=.false.
+! print *,"EVDW KURW",evdw,nres
+
+
+ return
+ end subroutine ecats_prot_amber
+
+!---------------------------------------------------------------------------
+! old for Ca2+
+ subroutine ecat_prot(ecation_prot)
+! use calc_data
+! use comm_momo
+ integer i,j,k,subchap,itmp,inum
+ real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,&
+ r7,r4,ecationcation
+ real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
+ dist_init,dist_temp,ecation_prot,rcal,rocal, &
+ Evan1,Evan2,EC,cm1mag,DASGL,delta,r0p,Epepcat, &
+ catl,cml,calpl, Etotal_p, Etotal_m,rtab,wdip,wmodquad,wquad1, &
+ wquad2,wvan1,E1,E2,wconst,wvan2,rcpm,dcmag,sin2thet,sinthet, &
+ costhet,v1m,v2m,wh2o,wc,rsecp,Ir,Irsecp,Irthrp,Irfourp,Irfiftp,&
+ Irsistp,Irseven,Irtwelv,Irthir,dE1dr,dE2dr,dEdcos,wquad2p,opt, &
+ rs,rthrp,rfourp,rsixp,reight,Irsixp,Ireight,Irtw,Irfourt, &
+ opt1,opt2,opt3,opt4,opt5,opt6,opt7,opt8,opt9,opt10,opt11,opt12,&
+ opt13,opt14,opt15,opt16,opt17,opt18,opt19, &
+ Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip,&
+ ndiv,ndivi
+ real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,&
+ gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, &
+ dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, &
+ tab1,tab2,tab3,diff,cm1,sc,p,tcat,talp,cm,drcp,drcp_norm,vcat, &
+ v1,v2,v3,myd_norm,dx,vcm,valpha,drdpep,dcosdpep,dcosddci,dEdpep,&
+ dEcCat,dEdipCm,dEdipCalp,dEquad1Cat,dEquad1Cm,dEquad1Calp, &
+ dEquad2Cat,dEquad2Cm,dEquad2Calpd,Evan1Cat,dEvan1Calp,dEvan2Cat,&
+ dEvan2Calp,dEtotalCat,dscvec,dEcCm,dEcCalp,dEdipCat,dEquad2Calp,&
+ dEvan1Cat
+ real(kind=8),dimension(6) :: vcatprm
+ ecation_prot=0.0d0
+! first lets calculate interaction with peptide groups
+ if (nres_molec(5).eq.0) return
+ itmp=0
+ do i=1,4
+ itmp=itmp+nres_molec(i)
+ enddo
+! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization
+ do i=ibond_start,ibond_end
+! cycle
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle ! leave dummy atoms
+ xi=0.5d0*(c(1,i)+c(1,i+1))
+ yi=0.5d0*(c(2,i)+c(2,i+1))
+ zi=0.5d0*(c(3,i)+c(3,i+1))
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+
+ do j=itmp+1,itmp+nres_molec(5)
+! print *,"WTF",itmp,j,i
+! all parameters were for Ca2+ to approximate single charge divide by two
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+ wconst=78*ndiv
+ wdip =1.092777950857032D2
+ wdip=wdip/wconst
+ wmodquad=-2.174122713004870D4
+ wmodquad=wmodquad/wconst
+ wquad1 = 3.901232068562804D1
+ wquad1=wquad1/wconst
+ wquad2 = 3
+ wquad2=wquad2/wconst
+ wvan1 = 0.1
+ wvan2 = 6
+! itmp=0
+
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+! enddo
+! enddo
+ rcpm = sqrt(xj**2+yj**2+zj**2)
+ drcp_norm(1)=xj/rcpm
+ drcp_norm(2)=yj/rcpm
+ drcp_norm(3)=zj/rcpm
+ dcmag=0.0
+ do k=1,3
+ dcmag=dcmag+dc(k,i)**2
+ enddo
+ dcmag=dsqrt(dcmag)
+ do k=1,3
+ myd_norm(k)=dc(k,i)/dcmag
+ enddo
+ costhet=drcp_norm(1)*myd_norm(1)+drcp_norm(2)*myd_norm(2)+&
+ drcp_norm(3)*myd_norm(3)
+ rsecp = rcpm**2
+ Ir = 1.0d0/rcpm
+ Irsecp = 1.0d0/rsecp
+ Irthrp = Irsecp/rcpm
+ Irfourp = Irthrp/rcpm
+ Irfiftp = Irfourp/rcpm
+ Irsistp=Irfiftp/rcpm
+ Irseven=Irsistp/rcpm
+ Irtwelv=Irsistp*Irsistp
+ Irthir=Irtwelv/rcpm
sin2thet = (1-costhet*costhet)
sinthet=sqrt(sin2thet)
E1 = wdip*Irsecp*costhet+(wmodquad*Irfourp+wquad1*Irthrp)&
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- &
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)
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
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
! 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)
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
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
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)
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)
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
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)
real(kind=8) :: facd4, adler, Fgb, facd3
integer troll,jj,istate
real (kind=8) :: dcosom1(3),dcosom2(3)
+ evdw=0.0d0
eps_out=80.0d0
sss_ele_cut=1.0d0
! print *,"EVDW KURW",evdw,nres
!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
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
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
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
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)
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!-------------------------------------------------------------------
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))
+ gradpepcatx(k,i) = gradpepcatx(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))
+! gradpepcatx(k,j) = gradpepcatx(k,j)+ dGCLdR * pom&
+! + dGGBdR * pom+ dGCVdR * pom&
+! + dPOLdR1 * (erhead_tail(k,1)&
+! -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j)))&
+! + dPOLdR2 * condor + dGLJdR * pom
+
+ gradpepcat(k,i) = gradpepcat(k,i) &
+ - dGCLdR * erhead(k)&
+ - dGGBdR * erhead(k)&
+ - dGCVdR * erhead(k)&
+ - dPOLdR1 * erhead_tail(k,1)&
+ - dPOLdR2 * erhead_tail(k,2)&
+ - dGLJdR * erhead(k)
+
+ gradpepcat(k,j) = gradpepcat(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
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
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) &
+ gradpepcatx(k,i) = gradpepcatx(k,i) &
- dPOLdR2 * (erhead_tail(k,2) &
-facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))
- gvdwx(k,j) = gvdwx(k,j) &
- + dPOLdR2 * condor
+! gradpepcatx(k,j) = gradpepcatx(k,j) &
+! + dPOLdR2 * condor
- gvdwc(k,i) = gvdwc(k,i) &
+ gradpepcat(k,i) = gradpepcat(k,i) &
- dPOLdR2 * erhead_tail(k,2)
- gvdwc(k,j) = gvdwc(k,j) &
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ dPOLdR2 * erhead_tail(k,2)
END DO
RETURN
- END SUBROUTINE enq
+ END SUBROUTINE enq_cat
+
SUBROUTINE eqd(Ecl,Elj,Epol)
use calc_data
use comm_momo
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)))
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx(k,i) = gvdwx(k,i) &
+ - dGCLdR * pom&
+ - dPOLdR1 * hawk &
+ - dGLJdR * pom
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx(k,j) = gvdwx(k,j) &
+ + dGCLdR * pom &
+ + dPOLdR1 * (erhead_tail(k,1) &
+ -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) &
+ + dGLJdR * pom
+
+
+ gvdwc(k,i) = gvdwc(k,i) &
+ - dGCLdR * erhead(k) &
+ - dPOLdR1 * erhead_tail(k,1) &
+ - dGLJdR * erhead(k)
+
+ gvdwc(k,j) = gvdwc(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1) &
+ + dGLJdR * erhead(k)
+
+ END DO
+ RETURN
+ END SUBROUTINE eqd
+ SUBROUTINE edq(Ecl,Elj,Epol)
+! IMPLICIT NONE
+ use comm_momo
+ use calc_data
+
+ double precision facd3, adler,ecl,elj,epol
+ alphapol2 = alphapol(itypj,itypi)
+ w1 = wqdip(1,itypi,itypj)
+ w2 = wqdip(2,itypi,itypj)
+ pis = sig0head(itypi,itypj)
+ eps_head = epshead(itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R2 = dsqrt(R2)
+
+!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj)
+!c! & +dhead(1,1,itypi,itypj))**2))
+!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj)
+!c! & +dhead(2,1,itypi,itypj))**2))
+
+
+!c!-------------------------------------------------------------------
+!c! ecl
+ sparrow = w1 * Qj * om1
+ hawk = w2 * Qj * Qj * (1.0d0 - sqom2)
+ ECL = sparrow / Rhead**2.0d0 &
+ - hawk / Rhead**4.0d0
+!c!-------------------------------------------------------------------
+!c! derivative of ecl is Gcl
+!c! dF/dr part
+ dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 &
+ + 4.0d0 * hawk / Rhead**5.0d0
+!c! dF/dom1
+ dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0)
+!c--------------------------------------------------------------------
+!c 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+nres) )
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
+ facd1 = d1 * vbld_inv(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)))
+
+ 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+nres))
+ 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
+
+ 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
+ write(iout,*) "KURWA2",Rhead
+ sparrow = w1 * Qj * om1
+ hawk = w2 * Qj * Qj * (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 * Qj) / (Rhead**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0)
+!c--------------------------------------------------------------------
+!c--------------------------------------------------------------------
+!c Polarization energy
+!c Epol
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR2 = R2 * R2 / MomoFac2
+ ee2 = exp(-(RR2 / (4.0d0 * a12sq)))
+ fgb2 = sqrt(RR2 + a12sq * ee2)
+ epol = 332.0d0 * eps_inout_fac * ((alphapol2/fgb2) ** 4.0d0 )
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0) &
+ / (fgb2 ** 5.0d0)
+ dFGBdR2 = ( (R2 / MomoFac2) &
+ * ( 2.0d0 - (0.5d0 * ee2) ) ) &
+ / (2.0d0 * fgb2)
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) &
+ * (2.0d0 - 0.5d0 * ee2) ) &
+ / (2.0d0 * fgb2)
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Elj
+ pom = (pis / Rhead)**6.0d0
+ Elj = 4.0d0 * eps_head * pom * (pom-1.0d0)
+!c! derivative of Elj is Glj
+ dGLJdR = 4.0d0 * eps_head &
+ * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) &
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))
+!c!-------------------------------------------------------------------
+
+!c! Return the results
+!c! (see comments in Eqq)
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ END DO
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j) )
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
+ facd1 = d1 * vbld_inv(i+nres)
+ facd2 = d2 * vbld_inv(j)
+ facd3 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres)
+ DO k = 1, 3
+ condor = (erhead_tail(k,2) &
+ + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j)))
pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
- gvdwx(k,i) = gvdwx(k,i) &
- - dGCLdR * pom&
- - dPOLdR1 * hawk &
- - dGLJdR * pom
+ gradpepcatx(k,i) = gradpepcatx(k,i) &
+ - dGCLdR * pom &
+ - dPOLdR2 * (erhead_tail(k,2) &
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) &
+ - dGLJdR * pom
- pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
- gvdwx(k,j) = gvdwx(k,j) &
- + dGCLdR * pom &
- + dPOLdR1 * (erhead_tail(k,1) &
- -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) &
- + dGLJdR * pom
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
+! gradpepcatx(k,j) = gradpepcatx(k,j) &
+! + dGCLdR * pom &
+! + dPOLdR2 * condor &
+! + dGLJdR * pom
- gvdwc(k,i) = gvdwc(k,i) &
- - dGCLdR * erhead(k) &
- - dPOLdR1 * erhead_tail(k,1) &
- - dGLJdR * erhead(k)
+ gradpepcat(k,i) = gradpepcat(k,i) &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k)
- gvdwc(k,j) = gvdwc(k,j) &
- + dGCLdR * erhead(k) &
- + dPOLdR1 * erhead_tail(k,1) &
- + dGLJdR * erhead(k)
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR2 * erhead_tail(k,2) &
+ + dGLJdR * erhead(k)
END DO
RETURN
- END SUBROUTINE eqd
- SUBROUTINE edq(Ecl,Elj,Epol)
-! IMPLICIT NONE
- use comm_momo
+ END SUBROUTINE edq_cat
+
+ SUBROUTINE edq_cat_pep(Ecl,Elj,Epol)
+ use comm_momo
use calc_data
double precision facd3, adler,ecl,elj,epol
- alphapol2 = alphapol(itypj,itypi)
- w1 = wqdip(1,itypi,itypj)
- w2 = wqdip(2,itypi,itypj)
- pis = sig0head(itypi,itypj)
- eps_head = epshead(itypi,itypj)
+ 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
!c!-------------------------------------------------------------------
!c! ecl
- sparrow = w1 * Qi * om1
- hawk = w2 * Qi * Qi * (1.0d0 - sqom2)
+ sparrow = w1 * Qj * om1
+ hawk = w2 * Qj * Qj * (1.0d0 - sqom2)
+! print *,"CO2", itypi,itypj
+! print *,"CO?!.", w1,w2,Qj,om1
ECL = sparrow / Rhead**2.0d0 &
- hawk / Rhead**4.0d0
!c!-------------------------------------------------------------------
dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 &
+ 4.0d0 * hawk / Rhead**5.0d0
!c! dF/dom1
- dGCLdOM1 = (w1 * Qi) / (Rhead**2.0d0)
+ dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0)
!c! dF/dom2
- dGCLdOM2 = (2.0d0 * w2 * Qi * Qi * om2) / (Rhead ** 4.0d0)
+ dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0)
+!c--------------------------------------------------------------------
!c--------------------------------------------------------------------
!c Polarization energy
!c Epol
* (((-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+nres) )
- eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) )
- adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
- facd1 = d1 * vbld_inv(i+nres)
- facd2 = d2 * vbld_inv(j+nres)
- facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+ erdxi = scalar( erhead(1), dC_norm(1,i) )
+ 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) )
+ facd1 = d1 * vbld_inv(i+1)/2.0
+ facd2 = d2 * vbld_inv(j)
+ facd3 = dtailcat(1,itypi,itypj) * vbld_inv(i+1)/2.0
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)))
- 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)+facd1*(erhead(k)-erdxi*dC_norm(k,i))
+! gradpepcatx(k,i) = gradpepcatx(k,i) &
+! - dGCLdR * pom &
+! - dPOLdR2 * (erhead_tail(k,2) &
+! -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) &
+! - dGLJdR * pom
- pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
- gvdwx(k,j) = gvdwx(k,j) &
- + dGCLdR * pom &
- + dPOLdR2 * condor &
- + dGLJdR * pom
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
+! gradpepcatx(k,j) = gradpepcatx(k,j) &
+! + dGCLdR * pom &
+! + dPOLdR2 * condor &
+! + dGLJdR * pom
- gvdwc(k,i) = gvdwc(k,i) &
+ gradpepcat(k,i) = gradpepcat(k,i) +0.5d0*( &
- dGCLdR * erhead(k) &
- dPOLdR2 * erhead_tail(k,2) &
- - dGLJdR * erhead(k)
+ - dGLJdR * erhead(k))
+ gradpepcat(k,i+1) = gradpepcat(k,i+1) +0.5d0*( &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k))
- gvdwc(k,j) = gvdwc(k,j) &
+
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ dGCLdR * erhead(k) &
+ dPOLdR2 * erhead_tail(k,2) &
+ dGLJdR * erhead(k)
END DO
RETURN
- END SUBROUTINE edq
+ END SUBROUTINE edq_cat_pep
+
SUBROUTINE edd(ECL)
! IMPLICIT NONE
use comm_momo
d1 = dhead(1, 1, itypi, itypj)
d2 = dhead(2, 1, itypi, itypj)
!c! ai*aj from Fgb
- a12sq = rborn(itypi,itypj) * rborn(itypj,itypi)
+ a12sq = rborn(itypi,itypj) * rborn(itypj,itypi)
+!c! a12sq = a12sq * a12sq
+!c! charge of amino acid itypi is...
+ Qi = icharge(itypi)
+ Qj = icharge(itypj)
+ Qij = Qi * Qj
+!c! chis1,2,12
+ chis1 = chis(itypi,itypj)
+ chis2 = chis(itypj,itypi)
+ chis12 = chis1 * chis2
+ sig1 = sigmap1(itypi,itypj)
+ sig2 = sigmap2(itypi,itypj)
+!c! write (*,*) "sig1 = ", sig1
+!c! write (*,*) "sig2 = ", sig2
+!c! alpha factors from Fcav/Gcav
+ b1cav = alphasur(1,itypi,itypj)
+! b1cav=0.0
+ b2cav = alphasur(2,itypi,itypj)
+ b3cav = alphasur(3,itypi,itypj)
+ b4cav = alphasur(4,itypi,itypj)
+ wqd = wquad(itypi, itypj)
+!c! used by Fgb
+ eps_in = epsintab(itypi,itypj)
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+!c! write (*,*) "eps_inout_fac = ", eps_inout_fac
+!c!-------------------------------------------------------------------
+!c! tail location and distance calculations
+ Rtail = 0.0d0
+ DO k = 1, 3
+ ctail(k,1)=c(k,i+nres)-dtail(1,itypi,itypj)*dc_norm(k,nres+i)
+ ctail(k,2)=c(k,j+nres)-dtail(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 = dhead(1, 1, itypi, itypj)
+ d2 = dhead(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
+
+
+ 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,5)
+!c! 1/(Gas Constant * Thermostate temperature) = BetaT
+!c! ENABLE THIS LINE WHEN USING CHECKGRAD!!!
+!c! t_bath = 300
+!c! BetaT = 1.0d0 / (t_bath * Rb)i
+ Rb=0.001986d0
+ BetaT = 1.0d0 / (298.0d0 * Rb)
+!c! Gay-berne var's
+ sig0ij = sigmacat( itypi,itypj )
+ chi1 = chi1cat( itypi, itypj )
+ chi2 = 0.0d0
+ chi12 = 0.0d0
+ chip1 = chipp1cat( itypi, itypj )
+ chip2 = 0.0d0
+ chip12 = 0.0d0
+!c! not used by momo potential, but needed by sc_angular which is shared
+!c! by all energy_potential subroutines
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+ dxj = 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 = rborn1cat(itypi,itypj) * rborn2cat(itypi,itypj)
+!c! a12sq = a12sq * a12sq
+!c! charge of amino acid itypi is...
+ Qi = icharge(itypi)
+ Qj = ichargecat(itypj)
+ Qij = Qi * Qj
+!c! chis1,2,12
+ chis1 = chis1cat(itypi,itypj)
+ chis2 = 0.0d0
+ chis12 = 0.0d0
+ sig1 = sigmap1cat(itypi,itypj)
+ sig2 = sigmap2cat(itypi,itypj)
+!c! alpha factors from Fcav/Gcav
+ b1cav = alphasurcat(1,itypi,itypj)
+ b2cav = alphasurcat(2,itypi,itypj)
+ b3cav = alphasurcat(3,itypi,itypj)
+ b4cav = alphasurcat(4,itypi,itypj)
+ wqd = wquadcat(itypi, itypj)
+!c! used by Fgb
+ eps_in = epsintabcat(itypi,itypj)
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+!c!-------------------------------------------------------------------
+!c! tail location and distance calculations
+ Rtail = 0.0d0
+ DO k = 1, 3
+ ctail(k,1)=c(k,i+nres)-dtailcat(1,itypi,itypj)*dc_norm(k,nres+i)
+ ctail(k,2)=c(k,j)!-dtailcat(2,itypi,itypj)*dc_norm(k,nres+j)
+ END DO
+!c! tail distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+ Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 )
+ Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 )
+ Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 )
+ Rtail = dsqrt( &
+ (Rtail_distance(1)*Rtail_distance(1)) &
+ + (Rtail_distance(2)*Rtail_distance(2)) &
+ + (Rtail_distance(3)*Rtail_distance(3)))
+!c!-------------------------------------------------------------------
+!c! Calculate location and distance between polar heads
+!c! distance between heads
+!c! for each one of our three dimensional space...
+ d1 = dheadcat(1, 1, itypi, itypj)
+ d2 = dheadcat(2, 1, itypi, itypj)
+
+ DO k = 1,3
+!c! location of polar head is computed by taking hydrophobic centre
+!c! and moving by a d1 * dc_norm vector
+!c! see unres publications for very informative images
+ chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres)
+ chead(k,2) = c(k, j)
+!c! distance
+!c! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+!c! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ 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
+
+ SUBROUTINE elgrad_init_cat_pep(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+ use comm_momo
+ use calc_data
+ real(kind=8) :: eheadtail,Egb,Ecl,Elj,Equad,Epol,Rb
+ eps_out=80.0d0
+ itypi = 10
+ itypj = itype(j,5)
+!c! 1/(Gas Constant * Thermostate temperature) = BetaT
+!c! ENABLE THIS LINE WHEN USING CHECKGRAD!!!
+!c! t_bath = 300
+!c! BetaT = 1.0d0 / (t_bath * Rb)i
+ Rb=0.001986d0
+ BetaT = 1.0d0 / (298.0d0 * Rb)
+!c! Gay-berne var's
+ sig0ij = sigmacat( itypi,itypj )
+ chi1 = chi1cat( itypi, itypj )
+ chi2 = 0.0d0
+ chi12 = 0.0d0
+ chip1 = chipp1cat( itypi, itypj )
+ chip2 = 0.0d0
+ chip12 = 0.0d0
+!c! not used by momo potential, but needed by sc_angular which is shared
+!c! by all energy_potential subroutines
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+ dxj = 0.0d0 !dc_norm( 1, nres+j )
+ dyj = 0.0d0 !dc_norm( 2, nres+j )
+ dzj = 0.0d0 !dc_norm( 3, nres+j )
+!c! distance from center of chain(?) to polar/charged head
+ d1 = dheadcat(1, 1, itypi, itypj)
+ d2 = dheadcat(2, 1, itypi, itypj)
+!c! ai*aj from Fgb
+ a12sq = rborn1cat(itypi,itypj) * rborn2cat(itypi,itypj)
!c! a12sq = a12sq * a12sq
!c! charge of amino acid itypi is...
- Qi = icharge(itypi)
- Qj = icharge(itypj)
- Qij = Qi * Qj
+ Qi = 0
+ Qj = ichargecat(itypj)
+! Qij = Qi * Qj
!c! chis1,2,12
- chis1 = chis(itypi,itypj)
- chis2 = chis(itypj,itypi)
- chis12 = chis1 * chis2
- sig1 = sigmap1(itypi,itypj)
- sig2 = sigmap2(itypi,itypj)
-!c! write (*,*) "sig1 = ", sig1
-!c! write (*,*) "sig2 = ", sig2
+ chis1 = chis1cat(itypi,itypj)
+ chis2 = 0.0d0
+ chis12 = 0.0d0
+ sig1 = sigmap1cat(itypi,itypj)
+ sig2 = sigmap2cat(itypi,itypj)
!c! alpha factors from Fcav/Gcav
- b1cav = alphasur(1,itypi,itypj)
-! b1cav=0.0
- b2cav = alphasur(2,itypi,itypj)
- b3cav = alphasur(3,itypi,itypj)
- b4cav = alphasur(4,itypi,itypj)
- wqd = wquad(itypi, itypj)
+ b1cav = alphasurcat(1,itypi,itypj)
+ b2cav = alphasurcat(2,itypi,itypj)
+ b3cav = alphasurcat(3,itypi,itypj)
+ b4cav = alphasurcat(4,itypi,itypj)
+ wqd = wquadcat(itypi, itypj)
!c! used by Fgb
- eps_in = epsintab(itypi,itypj)
+ eps_in = epsintabcat(itypi,itypj)
eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
-!c! write (*,*) "eps_inout_fac = ", eps_inout_fac
!c!-------------------------------------------------------------------
!c! tail location and distance calculations
Rtail = 0.0d0
DO k = 1, 3
- ctail(k,1)=c(k,i+nres)-dtail(1,itypi,itypj)*dc_norm(k,nres+i)
- ctail(k,2)=c(k,j+nres)-dtail(2,itypi,itypj)*dc_norm(k,nres+j)
+ ctail(k,1)=(c(k,i)+c(k,i+1))/2.0-dtailcat(1,itypi,itypj)*dc_norm(k,i)
+ ctail(k,2)=c(k,j)!-dtailcat(2,itypi,itypj)*dc_norm(k,nres+j)
END DO
!c! tail distances will be themselves usefull elswhere
!c1 (in Gcav, for example)
!c! Calculate location and distance between polar heads
!c! distance between heads
!c! for each one of our three dimensional space...
- d1 = dhead(1, 1, itypi, itypj)
- d2 = dhead(2, 1, itypi, itypj)
+ 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)
+ chead(k,1) = (c(k, i)+c(k,i+1))/2.0 + d1 * dc_norm(k, i)
+ chead(k,2) = c(k, j)
!c! distance
!c! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
!c! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
dPOLdOM1 = 0.0d0
dPOLdOM2 = 0.0d0
RETURN
- END SUBROUTINE elgrad_init
+ END SUBROUTINE elgrad_init_cat_pep
+
+ 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
+
+ subroutine make_SCSC_inter_list
+ include 'mpif.h'
+ real*8 :: xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp
+ real*8 :: dist_init, dist_temp,r_buff_list
+ integer:: contlisti(200*nres),contlistj(200*nres)
+! integer :: newcontlisti(200*nres),newcontlistj(200*nres)
+ integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_sc,g_ilist_sc
+ integer displ(0:nprocs),i_ilist_sc(0:nprocs),ierr
+! print *,"START make_SC"
+ r_buff_list=5.0
+ ilist_sc=0
+ do i=iatsc_s,iatsc_e
+ itypi=iabs(itype(i,1))
+ if (itypi.eq.ntyp1) cycle
+ 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=dmod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=dmod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ itypj=iabs(itype(j,1))
+ if (itypj.eq.ntyp1) cycle
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+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
+! r_buff_list is a read value for a buffer
+ if (sqrt(dist_init).le.(r_cut_ele+r_buff_list)) then
+! Here the list is created
+ ilist_sc=ilist_sc+1
+! this can be substituted by cantor and anti-cantor
+ contlisti(ilist_sc)=i
+ contlistj(ilist_sc)=j
+
+ endif
+ enddo
+ enddo
+ enddo
+! call MPI_Reduce(ilist_sc,g_ilist_sc,1,&
+! MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+! call MPI_Gather(newnss,1,MPI_INTEGER,&
+! i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
+#ifdef DEBUG
+ write (iout,*) "before MPIREDUCE",ilist_sc
+ do i=1,ilist_sc
+ write (iout,*) i,contlisti(i),contlistj(i)
+ enddo
+#endif
+ if (nfgtasks.gt.1)then
+
+ call MPI_Reduce(ilist_sc,g_ilist_sc,1,&
+ MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+ call MPI_Gather(ilist_sc,1,MPI_INTEGER,&
+ i_ilist_sc,1,MPI_INTEGER,king,FG_COMM,IERR)
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_ilist_sc(i-1)+displ(i-1)
+ enddo
+! write(iout,*) "before gather",displ(0),displ(1)
+ call MPI_Gatherv(contlisti,ilist_sc,MPI_INTEGER,&
+ newcontlisti,i_ilist_sc,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Gatherv(contlistj,ilist_sc,MPI_INTEGER,&
+ newcontlistj,i_ilist_sc,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+ call MPI_Bcast(newcontlisti,g_ilist_sc,MPI_INT,king,FG_COMM,IERR)
+ call MPI_Bcast(newcontlistj,g_ilist_sc,MPI_INT,king,FG_COMM,IERR)
+
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+
+ else
+ g_ilist_sc=ilist_sc
+
+ do i=1,ilist_sc
+ newcontlisti(i)=contlisti(i)
+ newcontlistj(i)=contlistj(i)
+ enddo
+ endif
+
+#ifdef DEBUG
+ write (iout,*) "after MPIREDUCE",g_ilist_sc
+ do i=1,g_ilist_sc
+ write (iout,*) i,newcontlisti(i),newcontlistj(i)
+ enddo
+#endif
+ call int_bounds(g_ilist_sc,g_listscsc_start,g_listscsc_end)
+ return
+ end subroutine make_SCSC_inter_list
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ subroutine make_SCp_inter_list
+ use MD_data, only: itime_mat
+
+ include 'mpif.h'
+ real*8 :: xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp
+ real*8 :: dist_init, dist_temp,r_buff_list
+ integer:: contlistscpi(200*nres),contlistscpj(200*nres)
+! integer :: newcontlistscpi(200*nres),newcontlistscpj(200*nres)
+ integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_scp,g_ilist_scp
+ integer displ(0:nprocs),i_ilist_scp(0:nprocs),ierr
+! print *,"START make_SC"
+ r_buff_list=5.0
+ ilist_scp=0
+ do i=iatscp_s,iatscp_e
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
+ 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 iint=1,nscp_gr(i)
+
+ do j=iscpstart(i,iint),iscpend(i,iint)
+ itypj=iabs(itype(j,1))
+ if (itypj.eq.ntyp1) cycle
+! Uncomment following three lines for SC-p interactions
+! xj=c(1,nres+j)-xi
+! yj=c(2,nres+j)-yi
+! zj=c(3,nres+j)-zi
+! Uncomment following three lines for Ca-p interactions
+! xj=c(1,j)-xi
+! yj=c(2,j)-yi
+! zj=c(3,j)-zi
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=mod(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
+#ifdef DEBUG
+ ! r_buff_list is a read value for a buffer
+ if ((sqrt(dist_init).le.(r_cut_ele)).and.(ifirstrun.eq.0)) then
+! Here the list is created
+ ilist_scp_first=ilist_scp_first+1
+! this can be substituted by cantor and anti-cantor
+ contlistscpi_f(ilist_scp_first)=i
+ contlistscpj_f(ilist_scp_first)=j
+ endif
+#endif
+! r_buff_list is a read value for a buffer
+ if (sqrt(dist_init).le.(r_cut_ele+r_buff_list)) then
+! Here the list is created
+ ilist_scp=ilist_scp+1
+! this can be substituted by cantor and anti-cantor
+ contlistscpi(ilist_scp)=i
+ contlistscpj(ilist_scp)=j
+ endif
+ enddo
+ enddo
+ enddo
+#ifdef DEBUG
+ write (iout,*) "before MPIREDUCE",ilist_scp
+ do i=1,ilist_scp
+ write (iout,*) i,contlistscpi(i),contlistscpj(i)
+ enddo
+#endif
+ if (nfgtasks.gt.1)then
+
+ call MPI_Reduce(ilist_scp,g_ilist_scp,1,&
+ MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+ call MPI_Gather(ilist_scp,1,MPI_INTEGER,&
+ i_ilist_scp,1,MPI_INTEGER,king,FG_COMM,IERR)
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_ilist_scp(i-1)+displ(i-1)
+ enddo
+! write(iout,*) "before gather",displ(0),displ(1)
+ call MPI_Gatherv(contlistscpi,ilist_scp,MPI_INTEGER,&
+ newcontlistscpi,i_ilist_scp,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Gatherv(contlistscpj,ilist_scp,MPI_INTEGER,&
+ newcontlistscpj,i_ilist_scp,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Bcast(g_ilist_scp,1,MPI_INT,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+ call MPI_Bcast(newcontlistscpi,g_ilist_scp,MPI_INT,king,FG_COMM,IERR)
+ call MPI_Bcast(newcontlistscpj,g_ilist_scp,MPI_INT,king,FG_COMM,IERR)
+
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+
+ else
+ g_ilist_scp=ilist_scp
+
+ do i=1,ilist_scp
+ newcontlistscpi(i)=contlistscpi(i)
+ newcontlistscpj(i)=contlistscpj(i)
+ enddo
+ endif
+
+#ifdef DEBUG
+ write (iout,*) "after MPIREDUCE",g_ilist_scp
+ do i=1,g_ilist_scp
+ write (iout,*) i,newcontlistscpi(i),newcontlistscpj(i)
+ enddo
+
+! if (ifirstrun.eq.0) ifirstrun=1
+! do i=1,ilist_scp_first
+! do j=1,g_ilist_scp
+! if ((newcontlistscpi(j).eq.contlistscpi_f(i)).and.&
+! (newcontlistscpj(j).eq.contlistscpj_f(i))) go to 126
+! enddo
+! print *,itime_mat,"ERROR matrix needs updating"
+! print *,contlistscpi_f(i),contlistscpj_f(i)
+! 126 continue
+! enddo
+#endif
+ call int_bounds(g_ilist_scp,g_listscp_start,g_listscp_end)
+
+ return
+ end subroutine make_SCp_inter_list
+
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+
+
+ subroutine make_pp_inter_list
+ include 'mpif.h'
+ real*8 :: xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp
+ real*8 :: xmedj,ymedj,zmedj
+ real*8 :: dist_init, dist_temp,r_buff_list,dxi,dyi,dzi,xmedi,ymedi,zmedi
+ real*8 :: dx_normi,dy_normi,dz_normi,dxj,dyj,dzj,dx_normj,dy_normj,dz_normj
+ integer:: contlistppi(200*nres),contlistppj(200*nres)
+! integer :: newcontlistppi(200*nres),newcontlistppj(200*nres)
+ integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_pp,g_ilist_pp
+ integer displ(0:nprocs),i_ilist_pp(0:nprocs),ierr
+! print *,"START make_SC"
+ ilist_pp=0
+ r_buff_list=5.0
+ do i=iatel_s,iatel_e
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
+ dxi=dc(1,i)
+ dyi=dc(2,i)
+ dzi=dc(3,i)
+ dx_normi=dc_norm(1,i)
+ dy_normi=dc_norm(2,i)
+ dz_normi=dc_norm(3,i)
+ xmedi=c(1,i)+0.5d0*dxi
+ ymedi=c(2,i)+0.5d0*dyi
+ zmedi=c(3,i)+0.5d0*dzi
+ xmedi=dmod(xmedi,boxxsize)
+ if (xmedi.lt.0) xmedi=xmedi+boxxsize
+ ymedi=dmod(ymedi,boxysize)
+ if (ymedi.lt.0) ymedi=ymedi+boxysize
+ zmedi=dmod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ do j=ielstart(i),ielend(i)
+! write (iout,*) i,j,itype(i,1),itype(j,1)
+ if (itype(j,1).eq.ntyp1.or. itype(j+1,1).eq.ntyp1) cycle
+
+! 1,j)
+ dxj=dc(1,j)
+ dyj=dc(2,j)
+ dzj=dc(3,j)
+ dx_normj=dc_norm(1,j)
+ dy_normj=dc_norm(2,j)
+ dz_normj=dc_norm(3,j)
+! xj=c(1,j)+0.5D0*dxj-xmedi
+! yj=c(2,j)+0.5D0*dyj-ymedi
+! zj=c(3,j)+0.5D0*dzj-zmedi
+ xj=c(1,j)+0.5D0*dxj
+ yj=c(2,j)+0.5D0*dyj
+ zj=c(3,j)+0.5D0*dzj
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+
+ dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ 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-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ endif
+ enddo
+ enddo
+ enddo
+
+ if (sqrt(dist_init).le.(r_cut_ele+r_buff_list)) then
+! Here the list is created
+ ilist_pp=ilist_pp+1
+! this can be substituted by cantor and anti-cantor
+ contlistppi(ilist_pp)=i
+ contlistppj(ilist_pp)=j
+ endif
+ enddo
+ enddo
+! enddo
+#ifdef DEBUG
+ write (iout,*) "before MPIREDUCE",ilist_pp
+ do i=1,ilist_pp
+ write (iout,*) i,contlistppi(i),contlistppj(i)
+ enddo
+#endif
+ if (nfgtasks.gt.1)then
+
+ call MPI_Reduce(ilist_pp,g_ilist_pp,1,&
+ MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+ call MPI_Gather(ilist_pp,1,MPI_INTEGER,&
+ i_ilist_pp,1,MPI_INTEGER,king,FG_COMM,IERR)
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_ilist_pp(i-1)+displ(i-1)
+ enddo
+! write(iout,*) "before gather",displ(0),displ(1)
+ call MPI_Gatherv(contlistppi,ilist_pp,MPI_INTEGER,&
+ newcontlistppi,i_ilist_pp,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Gatherv(contlistppj,ilist_pp,MPI_INTEGER,&
+ newcontlistppj,i_ilist_pp,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Bcast(g_ilist_pp,1,MPI_INT,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+ call MPI_Bcast(newcontlistppi,g_ilist_pp,MPI_INT,king,FG_COMM,IERR)
+ call MPI_Bcast(newcontlistppj,g_ilist_pp,MPI_INT,king,FG_COMM,IERR)
+
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+
+ else
+ g_ilist_pp=ilist_pp
+
+ do i=1,ilist_pp
+ newcontlistppi(i)=contlistppi(i)
+ newcontlistppj(i)=contlistppj(i)
+ enddo
+ endif
+ call int_bounds(g_ilist_pp,g_listpp_start,g_listpp_end)
+#ifdef DEBUG
+ write (iout,*) "after MPIREDUCE",g_ilist_pp
+ do i=1,g_ilist_pp
+ write (iout,*) i,newcontlistppi(i),newcontlistppj(i)
+ enddo
+#endif
+ return
+ end subroutine make_pp_inter_list
+
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+
+
+
end module energy
projects / unres4.git / blobdiff