From aedd2f72960d9cec69fcc7f0b6fd6555c7e22312 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Fri, 20 Nov 2015 10:19:56 +0100 Subject: [PATCH] adding of Czybyshev part 1 --- source/unres/src_MD-M/COMMON.CONTROL | 2 +- source/unres/src_MD-M/COMMON.DERIV | 2 +- source/unres/src_MD-M/COMMON.IOUNITS | 12 +- source/unres/src_MD-M/COMMON.TORSION | 14 +- source/unres/src_MD-M/DIMENSIONS | 3 +- source/unres/src_MD-M/energy_p_new-sep_barrier.F | 201 +++++++++++++++++++--- source/unres/src_MD-M/energy_p_new_barrier.F | 192 ++++++++++++++++++++- source/unres/src_MD-M/energy_split-sep.F | 50 +++++- source/unres/src_MD-M/initialize_p.F | 6 +- source/unres/src_MD-M/parmread.F | 27 +++ source/unres/src_MD-M/readrtns_CSA.F | 25 +++ 11 files changed, 498 insertions(+), 36 deletions(-) diff --git a/source/unres/src_MD-M/COMMON.CONTROL b/source/unres/src_MD-M/COMMON.CONTROL index 04cbfb7..a6cb608 100644 --- a/source/unres/src_MD-M/COMMON.CONTROL +++ b/source/unres/src_MD-M/COMMON.CONTROL @@ -1,6 +1,6 @@ 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, diff --git a/source/unres/src_MD-M/COMMON.DERIV b/source/unres/src_MD-M/COMMON.DERIV index 9b3b081..fcfe460 100644 --- a/source/unres/src_MD-M/COMMON.DERIV +++ b/source/unres/src_MD-M/COMMON.DERIV @@ -2,7 +2,7 @@ & 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), diff --git a/source/unres/src_MD-M/COMMON.IOUNITS b/source/unres/src_MD-M/COMMON.IOUNITS index 56e655e..8c4fdb7 100644 --- a/source/unres/src_MD-M/COMMON.IOUNITS +++ b/source/unres/src_MD-M/COMMON.IOUNITS @@ -11,11 +11,12 @@ C General I/O units & files 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 @@ -39,10 +40,13 @@ C CSA I/O units & files 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 diff --git a/source/unres/src_MD-M/COMMON.TORSION b/source/unres/src_MD-M/COMMON.TORSION index 4cf97e3..d9a2fea 100644 --- a/source/unres/src_MD-M/COMMON.TORSION +++ b/source/unres/src_MD-M/COMMON.TORSION @@ -1,15 +1,23 @@ 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 diff --git a/source/unres/src_MD-M/DIMENSIONS b/source/unres/src_MD-M/DIMENSIONS index fc17408..02deda8 100644 --- a/source/unres/src_MD-M/DIMENSIONS +++ b/source/unres/src_MD-M/DIMENSIONS @@ -48,8 +48,9 @@ C Number of AA types (at present only natural AA's will be handled 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, diff --git a/source/unres/src_MD-M/energy_p_new-sep_barrier.F b/source/unres/src_MD-M/energy_p_new-sep_barrier.F index c0b4c84..e0b1f44 100644 --- a/source/unres/src_MD-M/energy_p_new-sep_barrier.F +++ b/source/unres/src_MD-M/energy_p_new-sep_barrier.F @@ -1335,8 +1335,8 @@ C 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 @@ -1368,6 +1368,7 @@ C 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), @@ -1569,6 +1570,7 @@ C------------------------------------------------------------------------------- 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), @@ -1665,8 +1667,14 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions 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) @@ -1698,6 +1706,60 @@ C 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 @@ -1796,7 +1858,9 @@ C ggg(3)=facvdw*zj 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) @@ -1814,11 +1878,14 @@ cgrad enddo 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 @@ -2014,19 +2081,80 @@ cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij 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) @@ -2040,14 +2168,22 @@ cgrad enddo 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 @@ -2126,8 +2262,19 @@ c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2) 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 @@ -2195,17 +2342,29 @@ cgrad ghalfm=0.5D0*gggm(k) 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 diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index a889c75..8fdffe3 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -201,7 +201,14 @@ C 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 @@ -218,7 +225,14 @@ C Calculate the virtual-bond torsional energy. 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 @@ -227,7 +241,7 @@ c print *,"Processor",myrank," computed Utor" 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 @@ -7432,6 +7446,140 @@ C v1cij=v1c(1,j,itori,itori1,itori2,iblock,ntblock) 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 @@ -11171,4 +11319,46 @@ C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i) 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 diff --git a/source/unres/src_MD-M/energy_split-sep.F b/source/unres/src_MD-M/energy_split-sep.F index 24ab8dd..43f9baa 100644 --- a/source/unres/src_MD-M/energy_split-sep.F +++ b/source/unres/src_MD-M/energy_split-sep.F @@ -25,6 +25,7 @@ cMS$ATTRIBUTES C :: proc_proc 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 @@ -132,6 +133,10 @@ 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. @@ -261,6 +266,7 @@ cMS$ATTRIBUTES C :: proc_proc 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) @@ -393,6 +399,37 @@ C 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 @@ -414,7 +451,7 @@ C from other distance constraints. C C Calculate the virtual-bond-angle energy. C - call ebend(ebe) + call ebend(ebe,ethetcnstr) C C Calculate the SC local energy. C @@ -451,10 +488,19 @@ 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 diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F index da0d3f9..83ccf62 100644 --- a/source/unres/src_MD-M/initialize_p.F +++ b/source/unres/src_MD-M/initialize_p.F @@ -119,6 +119,8 @@ C 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 @@ -668,8 +670,8 @@ c & " ivec_start",ivec_start," ivec_end",ivec_end 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) diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index ae4d710..e156c97 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -775,6 +775,31 @@ C Martix of D parameters for two dimesional fourier series 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 @@ -1376,6 +1401,8 @@ C endif 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) diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index ebb59da..79840a4 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -155,6 +155,9 @@ C 2 reseved for further possible developement 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) @@ -2000,12 +2003,18 @@ C Get parameter filenames and open the parameter files. 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') @@ -2028,6 +2037,10 @@ c print *,"Processor",myrank," opened file IROTAM" 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') @@ -2040,6 +2053,8 @@ c print *,"Processor",myrank," opened file IFOURIER" 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) @@ -2057,6 +2072,10 @@ C Get parameter filenames and open the parameter files. 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) @@ -2065,6 +2084,8 @@ C Get parameter filenames and open the parameter files. 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) @@ -2081,6 +2102,10 @@ C Get parameter filenames and open the parameter files. 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 -- 1.7.9.5