integer modecalc,iscode,indpdb,indback,indphi,iranconf,icheckgrad,
& inprint,i2ndstr,mucadyn,constr_dist,symetr,AFMlog,selfguide,
- & shield_mode
+ & shield_mode,tor_mode
logical minim,refstr,pdbref,outpdb,outmol2,overlapsc,energy_dec,
& sideadd,lsecondary,read_cart,unres_pdb,
& vdisulf,searchsc,lmuca,dccart,extconf,out1file,
& gvdwpp,gel_loc,gel_loc_long,gvdwc_scpp,gliptranc,gliptranx,
& gradx_scp,gvdwc_scp,ghpbx,ghpbc,gloc,gloc_x,dtheta,dphi,dalpha,
& domega,gscloc,gsclocx,gradcorr,gradcorr_long,gradcorr5_long,
- & gradcorr6_long,gcorr6_turn_long,gvdwx,gshieldx
+ & gradcorr6_long,gcorr6_turn_long,gvdwx,gshieldx,gradafm
integer nfl,icg
common /derivat/ dcdv(6,maxdim),dxdv(6,maxdim),dxds(6,maxres),
& gradx(3,-1:maxres,2),gradc(3,-1:maxres,2),gvdwx(3,-1:maxres),
integer inp,iout,igeom,intin,ipdb,imol2,ipdbin,ithep,irotam,
& itorp,itordp,ifourier,ielep,isidep,iscpp,icbase,istat,
& ientin,ientout,izs1,isecpred,ibond,irest2,iifrag,icart,
- & irest1,isccor,ithep_pdb,irotam_pdb
+ & irest1,isccor,ithep_pdb,irotam_pdb,itorkcc,ithetkcc
common /iounits/ inp,iout,igeom,intin,ipdb,imol2,ipdbin,ithep,
& irotam,itorp,itordp,ifourier,ielep,isidep,iscpp,icbase,
& istat,ientin,ientout,izs1,isecpred,ibond,irest2,iifrag,
- & icart,irest1,isccor,ithep_pdb,irotam_pdb
+ & icart,irest1,isccor,ithep_pdb,irotam_pdb,
+ & itorkcc,ithetkcc
character*256 outname,intname,pdbname,mol2name,statname,intinname,
& entname,prefix,secpred,rest2name,qname,cartname,tmpdir,
& mremd_rst_name,curdir,pref_orig
C Parameter files
character*256 bondname,thetname,rotname,torname,tordname,
& fouriername,elename,sidename,scpname,sccorname,patname,
- & thetname_pdb,rotname_pdb,liptranname
+ & thetname_pdb,rotname_pdb,liptranname,thetkccname,
+ & torkccname
+
common /parfiles/ bondname,thetname,rotname,torname,tordname,
& fouriername,elename,sidename,scpname,sccorname,patname,
- & thetname_pdb,rotname_pdb,liptranname
+ & thetname_pdb,rotname_pdb,liptranname,thetkccname,
+ & torkccname
character*3 pot
C-----------------------------------------------------------------------
C INP - main input file
C Torsional constants of the rotation about virtual-bond dihedral angles
double precision v1,v2,vlor1,vlor2,vlor3,v0
- integer itortyp,ntortyp,nterm,nlor,nterm_old
+ integer itortyp,ntortyp,nterm,nlor,nterm_old,nterm_kcc_Tb,
+ & nterm_kcc,itortyp_kcc
common/torsion/v0(-maxtor:maxtor,-maxtor:maxtor,2),
& v1(maxterm,-maxtor:maxtor,-maxtor:maxtor,2),
& v2(maxterm,-maxtor:maxtor,-maxtor:maxtor,2),
& vlor1(maxlor,-maxtor:maxtor,-maxtor:maxtor),
& vlor2(maxlor,maxtor,maxtor),vlor3(maxlor,maxtor,maxtor),
+ & v1_kcc(maxtermkcc,-maxtor:maxtor,-maxtor:maxtor),
+ & v2_kcc(maxtermkcc,-maxtor:maxtor,-maxtor:maxtor),
+ & v1_chyb(maxtermkcc,-maxtor:maxtor,-maxtor:maxtor),
+ & v2_chyb(maxtermkcc,-maxtor:maxtor,-maxtor:maxtor),
& itortyp(-ntyp1:ntyp1),ntortyp,
+ & itortyp_kcc(-ntyp1:ntyp1),
& nterm(-maxtor:maxtor,-maxtor:maxtor,2),
- & nlor(-maxtor:maxtor,-maxtor:maxtor,2)
- & ,nterm_old
+ & nlor(-maxtor:maxtor,-maxtor:maxtor,2),
+ & nterm_kcc_Tb(-maxtor:maxtor,-maxtor:maxtor),
+ & nterm_kcc(-maxtor:maxtor,-maxtor:maxtor),
+ & nterm_old
C 6/23/01 - constants for double torsionals
double precision v1c,v1s,v2c,v2s
integer ntermd_1,ntermd_2
parameter (ntyp=24,ntyp1=ntyp+1)
C Max. number of types of dihedral angles & multiplicity of torsional barriers
C and the number of terms in double torsionals
- integer maxtor,maxterm,maxlor,maxtermd_1,maxtermd_2
+ integer maxtor,maxterm,maxlor,maxtermd_1,maxtermd_2,maxtermkcc
parameter (maxtor=4,maxterm=10,maxlor=3,maxtermd_1=8,maxtermd_2=8)
+ parameter (maxtermkcc=20)
C Max. number of residue types and parameters in expressions for
C virtual-bond angle bending potentials
integer maxthetyp,maxthetyp1,maxtheterm,maxtheterm2,maxtheterm3,
C Calculate the components of the gradient in DC and X
C
do l=1,3
- gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k)
- gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k)
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
enddo
return
end
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
include 'COMMON.TIME1'
+ include 'COMMON.SHIELD'
dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
& erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
include 'COMMON.TIME1'
+ include 'COMMON.SHIELD'
dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
& erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
fac3=ael6i*r6ij
fac4=ael3i*r3ij
evdwij=ev1+ev2
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+ endif
el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
el2=fac4*fac
+ el1=el1*fac_shield(i)**2*fac_shield(j)**2
+ el2=el2*fac_shield(i)**2*fac_shield(j)**2
eesij=el1+el2
C 12/26/95 - for the evaluation of multi-body H-bonding interactions
ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
ggg(1)=facel*xj
ggg(2)=facel*yj
ggg(3)=facel*zj
+ if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+ & (shield_mode.gt.0)) then
+C print *,i,j
+ do ilist=1,ishield_list(i)
+ iresshield=shield_list(ilist,i)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
+ & *2.0
+ gshieldx(k,iresshield)=gshieldx(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0
+ gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
+C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
+C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C if (iresshield.gt.i) then
+C do ishi=i+1,iresshield-1
+C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
+C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C
+C enddo
+C else
+C do ishi=iresshield,i
+C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
+C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C
+C enddo
+C endif
+ enddo
+ enddo
+ do ilist=1,ishield_list(j)
+ iresshield=shield_list(ilist,j)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
+ & *2.0
+ gshieldx(k,iresshield)=gshieldx(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
+ gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
+ enddo
+ enddo
+
+ do k=1,3
+ gshieldc(k,i)=gshieldc(k,i)+
+ & grad_shield(k,i)*eesij/fac_shield(i)*2.0
+ gshieldc(k,j)=gshieldc(k,j)+
+ & grad_shield(k,j)*eesij/fac_shield(j)*2.0
+ gshieldc(k,i-1)=gshieldc(k,i-1)+
+ & grad_shield(k,i)*eesij/fac_shield(i)*2.0
+ gshieldc(k,j-1)=gshieldc(k,j-1)+
+ & grad_shield(k,j)*eesij/fac_shield(j)*2.0
+
+ enddo
+ endif
+
c do k=1,3
c ghalf=0.5D0*ggg(k)
c gelc(k,i)=gelc(k,i)+ghalf
cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
cd & (dcosg(k),k=1,3)
do k=1,3
- ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
+ ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
+ & fac_shield(i)**2*fac_shield(j)**2
+
enddo
c do k=1,3
c ghalf=0.5D0*ggg(k)
cgrad enddo
do k=1,3
gelc(k,i)=gelc(k,i)
- & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
- & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+ & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
+ & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
+ & *fac_shield(i)**2*fac_shield(j)**2
+
gelc(k,j)=gelc(k,j)
- & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
- & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+ & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
+ & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
+ & *fac_shield(i)**2*fac_shield(j)**2
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
enddo
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
& 'eelloc',i,j,eel_loc_ij
+
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+C else
+C fac_shield(i)=0.4
+C fac_shield(j)=0.6
+ endif
+ eel_loc_ij=eel_loc_ij
+ & *fac_shield(i)*fac_shield(j)
eel_loc=eel_loc+eel_loc_ij
+
+ if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+ & (shield_mode.gt.0)) then
+C print *,i,j
+
+ do ilist=1,ishield_list(i)
+ iresshield=shield_list(ilist,i)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
+ & /fac_shield(i)
+C & *2.0
+ gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
+ gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
+ & +rlocshield
+ enddo
+ enddo
+ do ilist=1,ishield_list(j)
+ iresshield=shield_list(ilist,j)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
+ & /fac_shield(j)
+C & *2.0
+ gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
+ gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
+ & +rlocshield
+
+ enddo
+ enddo
+
+ do k=1,3
+ gshieldc_ll(k,i)=gshieldc_ll(k,i)+
+ & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
+ gshieldc_ll(k,j)=gshieldc_ll(k,j)+
+ & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
+ gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
+ & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
+ gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
+ & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
+ enddo
+ endif
+
C Partial derivatives in virtual-bond dihedral angles gamma
if (i.gt.1)
& gel_loc_loc(i-1)=gel_loc_loc(i-1)+
- & a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
- & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
+ & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
+ & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
+ & *fac_shield(i)*fac_shield(j)
+
gel_loc_loc(j-1)=gel_loc_loc(j-1)+
- & a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
- & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
+ & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
+ & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
+ & *fac_shield(i)*fac_shield(j)
+
C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
do l=1,3
- ggg(l)=agg(l,1)*muij(1)+
- & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
+ ggg(l)=(agg(l,1)*muij(1)+
+ & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
+ & *fac_shield(i)*fac_shield(j)
+
gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
cgrad ghalf=0.5d0*ggg(l)
cgrad enddo
C Remaining derivatives of eello
do l=1,3
- gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
- & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
- gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
- & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
- gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
- & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
- gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
- & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
+ gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
+ & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
+ & *fac_shield(i)*fac_shield(j)
+
+ gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
+ & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
+ & *fac_shield(i)*fac_shield(j)
+
+ gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
+ & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
+ & *fac_shield(i)*fac_shield(j)
+
+ gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
+ & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
+ & *fac_shield(i)*fac_shield(j)
+
enddo
ENDIF
C Change 12/26/95 to calculate four-body contributions to H-bonding energy
ees0mij=0
endif
c ees0mij=0.0D0
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0d0
+ fac_shield(j)=1.0d0
+ else
+ ees0plist(num_conti,i)=j
+C fac_shield(i)=0.4d0
+C fac_shield(j)=0.6d0
+ endif
ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
+ & *fac_shield(i)*fac_shield(j)
ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
+ & *fac_shield(i)*fac_shield(j)
+
C Diagnostics. Comment out or remove after debugging!
c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
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)
+ & *fac_shield(i)*fac_shield(j)
+
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)
+ & *fac_shield(i)*fac_shield(j)
+
gacontp_hb3(k,num_conti,i)=gggp(k)
+ & *fac_shield(i)*fac_shield(j)
+
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)
+ & *fac_shield(i)*fac_shield(j)
+
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)
+ & *fac_shield(i)*fac_shield(j)
+
gacontm_hb3(k,num_conti,i)=gggm(k)
+ & *fac_shield(i)*fac_shield(j)
+
enddo
ENDIF ! wcorr
endif ! num_conti.le.maxconts
C Calculate the virtual-bond-angle energy.
C
if (wang.gt.0d0) then
+ if ((tor_mode.eq.0).or.(tor_mode.eq.2)) then
call ebend(ebe,ethetacnstr)
+ endif
+C ebend kcc is Kubo cumulant clustered rigorous attemp to derive the
+C energy function
+ if ((tor_mode.eq.1).or.(tor_mode.eq.2)) then
+ call ebend_kcc(ebe,ethetacnstr)
+ endif
else
ebe=0
ethetacnstr=0
C
cd print *,'nterm=',nterm
if (wtor.gt.0) then
+ if ((tor_mode.eq.0).or.(tor_mode.eq.2)) then
call etor(etors,edihcnstr)
+ endif
+C etor kcc is Kubo cumulant clustered rigorous attemp to derive the
+C energy function
+ if ((tor_mode.eq.1).or.(tor_mode.eq.2)) then
+ call etor(etors,edihcnstr)
+ endif
else
etors=0
edihcnstr=0
C
C 6/23/01 Calculate double-torsional energy
C
- if (wtor_d.gt.0) then
+ if ((wtor_d.gt.0).and.((tor_mode.eq.0).or.(tor_mode.eq.2))) then
call etor_d(etors_d)
else
etors_d=0
return
end
#endif
+C----------------------------------------------------------------------------------
+C The rigorous attempt to derive energy function
+ subroutine etor_kcc(etors,edihcnstr)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.VAR'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.TORSION'
+ include 'COMMON.INTERACT'
+ include 'COMMON.DERIV'
+ include 'COMMON.CHAIN'
+ include 'COMMON.NAMES'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.FFIELD'
+ include 'COMMON.TORCNSTR'
+ include 'COMMON.CONTROL'
+ logical lprn
+ double precision thybt1(maxtermkcc),thybt2(maxtermkcc)
+C Set lprn=.true. for debugging
+ lprn=.false.
+c lprn=.true.
+ 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).eq.ntyp1.or. itype(i-1).eq.ntyp1
+ & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
+ itori=itortyp_kcc(itype(i-2))
+ itori1=itortyp_kcc(itype(i-1))
+ phii=phi(i)
+ glocig=0.0D0
+ glocit1=0.0d0
+ glocit2=0.0d0
+ sumnonchebyshev=0.0d0
+ sumchebyshev=0.0d0
+C to avoid multiple devision by 2
+ theti22=0.5d0*theta(i)
+C theta 12 is the theta_1 /2
+C theta 22 is theta_2 /2
+ theti12=0.5d0*theta(i-1)
+C and appropriate sinus function
+ sinthet2=dsin(theta(i))
+ sinthet1=dsin(theta(i-1))
+ costhet1=dcos(theta(i-1))
+ costhet2=dcos(theta(i))
+C to speed up lets store its mutliplication
+ sint1t2=sinthet2*sinthet1
+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
+ do j=1,nterm_kcc(itori,itori1)
+
+ v1ij=v1_kcc(j,itori,itori1)
+ v2ij=v2_kcc(j,itori,itori1)
+C v1ij is c_n and d_n in euation above
+ cosphi=dcos(j*phii)
+ sinphi=dsin(j*phii)
+ sint1t2n=sint1t2**j
+ sumnonchebyshev=sumnonchebyshev+
+ & sint1t2n*(v1ij*cosphi+v2ij*sinphi)
+ actval=sint1t2n*(v1ij*cosphi+v2ij*sinphi)
+C etors=etors+sint1t2n*(v1ij*cosphi+v2ij*sinphi)
+C if (energy_dec) etors_ii=etors_ii+
+C & v1ij*cosphi+v2ij*sinphi
+C glocig is the gradient local i site in gamma
+ glocig=glocig+j*(v2ij*cosphi-v1ij*sinphi)*sint1t2n
+C now gradient over theta_1
+ glocit1=glocit1+actval/sinthet1*j*costhet1
+ glocit2=glocit2+actval/sinthet2*j*costhet2
+ enddo
+
+C now the Czebyshev polinominal sum
+ do j=1,nterm_kcc_Tb(itori,itori1)
+ thybt1(j)=v1_chyb(j,itori,itori1)
+ thybt2(j)=v2_chyb(j,itori,itori1)
+ enddo
+ sumth1thyb=tschebyshev
+ & (1,nterm_kcc_Tb(itori,itori1),thybt1(1),dcos(theta12))
+ gradthybt1=gradtschebyshev
+ & (0,nterm_kcc_Tb(itori,itori1)-1,thybt1(1),dcos(theta12))
+ & *(nterm_kcc_Tb(itori,itori1))*0.5*dsin(theta12)
+ sumth2thyb=tschebyshev
+ & (1,nterm_kcc_Tb(itori,itori1),thybt2(1),dcos(theta22))
+ gradthybt2=gradtschebyshev
+ & (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1),dcos(theta22))
+ & *(nterm_kcc_Tb(itori,itori1))*0.5*dsin(theta22)
+
+C now overal sumation
+ etors=etors+(1.0d0+sumth1thyb+sumth2thyb)*sumnonchebyshev
+C derivative over gamma
+ gloc(i-3,icg)=gloc(i-3,icg)+wtor*glocig
+ & *(1.0d0+sumth1thyb+sumth2thyb)
+C derivative over theta1
+ gloc(nphi+i-3,icg)=gloc(nphi+i-3,icg)+wang*
+ & (glocit1*(1.0d0+sumth1thyb+sumth2thyb)+
+ & sumnonchebyshev*gradthybt1)
+C now derivative over theta2
+ gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*
+ & (glocit2*(1.0d0+sumth1thyb+sumth2thyb)+
+ & sumnonchebyshev*gradthybt2)
+ enddo
+
+
+C gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
+! 6/20/98 - dihedral angle constraints
+ if (tor_mode.ne.2) then
+ edihcnstr=0.0d0
+c do i=1,ndih_constr
+ do i=idihconstr_start,idihconstr_end
+ itori=idih_constr(i)
+ phii=phi(itori)
+ difi=pinorm(phii-phi0(i))
+ if (difi.gt.drange(i)) then
+ difi=difi-drange(i)
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
+ else if (difi.lt.-drange(i)) then
+ difi=difi+drange(i)
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
+ else
+ difi=0.0
+ endif
+ enddo
+ endif
+ return
+ end
+
+
+
c------------------------------------------------------------------------------
subroutine eback_sc_corr(esccor)
c 7/21/2007 Correlations between the backbone-local and side-chain-local
enddo
return
end
+C--------------------------------------------------------------------------
+ double precision function tschebyshev(m,n,x,y)
+ implicit none
+ include "DIMENSIONS"
+ 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
+C--------------------------------------------------------------------------
+ double precision function gradtschebyshev(m,n,x,y)
+ implicit none
+ include "DIMENSIONS"
+ 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)=2.0d0*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
+ gradtschebyshev=aux
+ return
+ end
include 'COMMON.VAR'
include 'COMMON.LOCAL'
include 'COMMON.MD'
+ include 'COMMON.CONTROL'
c write(iout,'(a,i2)')'Calling etotal_long ipot=',ipot
if (modecalc.eq.12.or.modecalc.eq.14) then
#ifdef MPI
C
107 continue
call vec_and_deriv
+ if (shield_mode.gt.0) then
+ call set_shield_fac
+ endif
+
if (ipot.lt.6) then
#ifdef SPLITELE
if (welec.gt.0d0.or.wvdwpp.gt.0d0.or.wel_loc.gt.0d0.or.
include 'COMMON.CHAIN'
include 'COMMON.VAR'
include 'COMMON.LOCAL'
+ include 'COMMON.CONTROL'
c write(iout,'(a,i2)')'Calling etotal_short ipot=',ipot
c call flush(iout)
C Calculate electrostatic (H-bonding) energy of the main chain.
C
107 continue
+ call vec_and_deriv
+ if (shield_mode.gt.0) then
+ call set_shield_fac
+ endif
+
+ if (ipot.lt.6) then
+#ifdef SPLITELE
+ if (welec.gt.0d0.or.wvdwpp.gt.0d0.or.wel_loc.gt.0d0.or.
+ & wturn3.gt.0d0.or.wturn4.gt.0d0 .or. wcorr.gt.0.0d0
+ & .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0
+ & .or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then
+#else
+ if (welec.gt.0d0.or.wel_loc.gt.0d0.or.
+ & wturn3.gt.0d0.or.wturn4.gt.0d0 .or. wcorr.gt.0.0d0
+ & .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0
+ & .or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then
+#endif
+ call eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
+ else
+ ees=0
+ evdw1=0
+ eel_loc=0
+ eello_turn3=0
+ eello_turn4=0
+ endif
+ else
+c write (iout,*) "Soft-spheer ELEC potential"
+ call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
+ & eello_turn4)
+ endif
+
c
c Calculate the short-range part of Evdwpp
c
C
C Calculate the virtual-bond-angle energy.
C
- call ebend(ebe)
+ call ebend(ebe,ethetcnstr)
C
C Calculate the SC local energy.
C
energia(18)=0.0d0
#endif
#ifdef SPLITELE
+ energia(3)=ees
energia(16)=evdw1
#else
- energia(3)=evdw1
+ energia(3)=ees+evdw1
+ energia(16)=0.0d0
#endif
+ energia(4)=ecorr
+ energia(5)=ecorr5
+ energia(6)=ecorr6
+ energia(7)=eel_loc
+ energia(8)=eello_turn3
+ energia(9)=eello_turn4
+
energia(11)=ebe
energia(12)=escloc
energia(13)=etors
icsa_in=40
crc for ifc error 118
icsa_pdb=42
+ itorkcc=43
+ ithetkcc=44
C Lipidic input file for parameters range 60-79
iliptranpar=60
C input file for transfer sidechain and peptide group inside the
call int_bounds(ndih_constr,idihconstr_start,idihconstr_end)
endif
if (ntheta_constr.eq.0) then
- idihconstr_start=1
- idihconstr_end=0
+ ithetaconstr_start=1
+ ithetaconstr_end=0
else
call int_bounds
& (ntheta_constr,ithetaconstr_start,ithetaconstr_end)
enddo
endif
#endif
+C read Czybyshev torsional parameters
+ read (itorkcc,*,end=121,err=121) nkcctyp
+ read (itorkcc,*,end=121,err=121) (itortyp_kcc(i),i=1,ntyp)
+ do i=-ntyp,-1
+ itortyp_kcc(i)=-itortyp_kcc(-i)
+ enddo
+ do i=1,nkcctyp
+ do j=1,nkcctyp
+C first we read the cos and sin gamma parameters
+ read (itorkcc,*,end=121,err=121) nterm_kcc(j,i)
+ do k=1,nterm_kcc(j,i)
+ read (itorkcc,*,end=121,err=121)
+ & v1_kcc(k,j,i),v2_kcc(k,j,i)
+ enddo
+C now Czybyshev parameters
+ read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
+ do k=1,nterm_kcc_Tb(j,i)
+ read (itorkcc,*,end=121,err=121)
+ & v1_chyb(k,j,i),v2_chyb(k,j,i)
+ enddo
+ enddo
+ enddo
+C here will be the apropriate recalibrating for D-aminoacid
+
+
C Read of Side-chain backbone correlation parameters
C Modified 11 May 2012 by Adasko
CCC
118 write (iout,*) "Error reading SCp interaction parameters."
goto 999
119 write (iout,*) "Error reading SCCOR parameters"
+ go to 999
+ 121 write (iout,*) "Error in Czybyshev parameters"
999 continue
#ifdef MPI
call MPI_Finalize(Ierror)
C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
write(iout,*) "shield_mode",shield_mode
C endif
+ call readi(controlcard,'TORMODE',tor_mode,0)
+C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+ write(iout,*) "torsional and valence angle mode",tor_mode
call readi(controlcard,'MAXGEN',maxgen,10000)
call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
call readi(controlcard,"KDIAG",kdiag,0)
open (itorp,file=torname,status='old',readonly,shared)
call getenv_loc('TORDPAR',tordname)
open (itordp,file=tordname,status='old',readonly,shared)
+ call getenv_loc('TORKCC',torkccname)
+ open (itorkcc,file=torkccname,status='old',readonly,shared)
+ call getenv_loc('THETKCC',thetkccname)
+ open (ithetkcc,file=thetkccname,status='old',readonly,shared)
call getenv_loc('FOURIER',fouriername)
open (ifourier,file=fouriername,status='old',readonly,shared)
call getenv_loc('ELEPAR',elename)
open (ielep,file=elename,status='old',readonly,shared)
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old',readonly,shared)
+ call getenv_loc('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old',readonly,shared)
#elif (defined CRAY) || (defined AIX)
open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
& action='read')
c print *,"Processor",myrank," opened file ITORP"
call getenv_loc('TORDPAR',tordname)
open (itordp,file=tordname,status='old',action='read')
+ call getenv_loc('TORKCC',torkccname)
+ open (itorkcc,file=torkccname,status='old',action='read')
+ call getenv_loc('THETKCC',thetkccname)
+ open (ithetkcc,file=thetkccname,status='old',action='read')
c print *,"Processor",myrank," opened file ITORDP"
call getenv_loc('SCCORPAR',sccorname)
open (isccor,file=sccorname,status='old',action='read')
c print *,"Processor",myrank," opened file IELEP"
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old',action='read')
+ call getenv_loc('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old',action='read')
c print *,"Processor",myrank," opened file ISIDEP"
c print *,"Processor",myrank," opened parameter files"
#elif (defined G77)
open (itorp,file=torname,status='old')
call getenv_loc('TORDPAR',tordname)
open (itordp,file=tordname,status='old')
+ call getenv_loc('TORKCC',torkccname)
+ open (itorkcc,file=torkccname,status='old')
+ call getenv_loc('THETKCC',thetkccname)
+ open (ithetkcc,file=thetkccname,status='old')
call getenv_loc('SCCORPAR',sccorname)
open (isccor,file=sccorname,status='old')
call getenv_loc('FOURIER',fouriername)
open (ielep,file=elename,status='old')
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old')
+ call getenv_loc('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old')
#else
open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
& readonly)
open (itorp,file=torname,status='old',readonly)
call getenv_loc('TORDPAR',tordname)
open (itordp,file=tordname,status='old',readonly)
+ call getenv_loc('TORKCC',torkccname)
+ open (itorkcc,file=torkccname,status='old',readonly)
+ call getenv_loc('THETKCC',thetkccname)
+ open (ithetkcc,file=thetkccname,status='old',readonly)
call getenv_loc('SCCORPAR',sccorname)
open (isccor,file=sccorname,status='old',readonly)
#ifndef CRYST_THETA
projects / unres.git / commitdiff