From 24d71016d39facb439708acae299529762f51e81 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Thu, 11 Feb 2016 07:46:01 +0100 Subject: [PATCH] introduction of shielding to cluster DEBUG mode --- source/cluster/wham/src-M/COMMON.CLUSTER | 8 +- source/cluster/wham/src-M/COMMON.CONTROL | 3 +- source/cluster/wham/src-M/COMMON.DERIV | 20 +- source/cluster/wham/src-M/DIMENSIONS | 4 +- source/cluster/wham/src-M/energy_p_new.F | 907 +++++++++++++++++++++++++- source/cluster/wham/src-M/main_clust.F | 4 +- source/cluster/wham/src-M/parmread.F | 10 +- source/cluster/wham/src-M/probabl.F | 22 +- source/cluster/wham/src-M/read_coords.F | 12 +- source/cluster/wham/src-M/readrtns.F | 30 +- source/cluster/wham/src-M/work_partition.F | 13 +- source/unres/src_MD-M/energy_p_new_barrier.F | 1 + source/wham/src-M/energy_p_new.F | 51 +- source/wham/src-M/parmread.F | 2 +- source/wham/src-M/wham_calc1.F | 8 +- 15 files changed, 1031 insertions(+), 64 deletions(-) diff --git a/source/cluster/wham/src-M/COMMON.CLUSTER b/source/cluster/wham/src-M/COMMON.CLUSTER index 4477d19..4fedc52 100644 --- a/source/cluster/wham/src-M/COMMON.CLUSTER +++ b/source/cluster/wham/src-M/COMMON.CLUSTER @@ -6,12 +6,12 @@ integer ncut,ngr,licz,nconf,iass,icc,mult,list_conf, & nss_all,ihpb_all,jhpb_all,iass_tot,iscore,nprop common /clu/ diss(maxdist),energy(0:maxconf), - & enetb(max_ene,maxstr_proc),ecut, + & enetb(max_ene,maxconf),ecut, & entfac(maxconf),totfree(0:maxconf),totfree_gr(maxgr), & rcutoff(max_cut+1),ncut,min_var,tree,plot_tree,lgrp common /clu1/ ngr,licz(maxgr),nconf(maxgr,maxingr),iass(maxgr), & iass_tot(maxgr,max_cut),list_conf(maxconf) - common /alles/ allcart(3,maxres2,maxstr_proc),rmstb(maxconf), + common /alles/ allcart(3,maxres2,maxconf),rmstb(maxconf), & icc(maxconf), - & mult(maxres),nss_all(maxstr_proc),ihpb_all(maxss,maxstr_proc), - & jhpb_all(maxss,maxstr_proc),iscore(maxconf),nprop + & mult(maxres),nss_all(maxconf),ihpb_all(maxss,maxconf), + & jhpb_all(maxss,maxconf),iscore(maxconf),nprop diff --git a/source/cluster/wham/src-M/COMMON.CONTROL b/source/cluster/wham/src-M/COMMON.CONTROL index f339ae9..708fe2b 100644 --- a/source/cluster/wham/src-M/COMMON.CONTROL +++ b/source/cluster/wham/src-M/COMMON.CONTROL @@ -1,6 +1,6 @@ double precision betaT integer iscode,indpdb,outpdb,outmol2,iopt,nstart,nend,symetr, - & constr_dist + & constr_dist,shield_mode,tor_mode logical refstr,pdbref,punch_dist,print_dist,caonly,lside, & lprint_cart,lprint_int,from_cart,efree,from_bx,from_cx, & with_dihed_constr,with_theta_constr @@ -8,4 +8,5 @@ & punch_dist,print_dist,caonly,lside,lprint_cart,lprint_int, & from_cart,from_bx,from_cx, with_dihed_constr,with_theta_constr, & efree,iopt,nstart,nend,symetr, + & tor_mode,shield_mode, & constr_dist diff --git a/source/cluster/wham/src-M/COMMON.DERIV b/source/cluster/wham/src-M/COMMON.DERIV index 79f8630..4e116e1 100644 --- a/source/cluster/wham/src-M/COMMON.DERIV +++ b/source/cluster/wham/src-M/COMMON.DERIV @@ -3,13 +3,31 @@ & gradcorr5,gradcorr6,gel_loc,gcorr3_turn,gcorr4_turn,gcorr6_turn, & gel_loc_loc,gel_loc_turn3,gel_loc_turn4,gel_loc_turn6,gcorr_loc, & g_corr5_loc,g_corr6_loc,gradb,gradbx,gsccorc,gsccorx,gsccor_loc, - & gscloc,gsclocx + & gscloc,gsclocx,gshieldx,gradafm, + & gshieldc, gshieldc_loc, gshieldx_ec, gshieldc_ec, + & gshieldc_loc_ec, gshieldx_t3,gshieldc_t3,gshieldc_loc_t3, + & gshieldx_t4, gshieldc_t4,gshieldc_loc_t4,gshieldx_ll, + & gshieldc_ll, gshieldc_loc_ll + + integer nfl,icg logical calc_grad common /derivat/ dcdv(6,maxdim),dxdv(6,maxdim),dxds(6,maxres), & gradx(3,maxres,2),gradc(3,maxres,2),gvdwx(3,maxres), & gvdwc(3,maxres),gelc(3,maxres),gvdwpp(3,maxres), & gradx_scp(3,maxres), + & gliptranc(3,-1:maxres), + & gliptranx(3,-1:maxres), + & gshieldx(3,-1:maxres), gshieldc(3,-1:maxres), + & gshieldc_loc(3,-1:maxres), + & gshieldx_ec(3,-1:maxres), gshieldc_ec(3,-1:maxres), + & gshieldc_loc_ec(3,-1:maxres), + & gshieldx_t3(3,-1:maxres), gshieldc_t3(3,-1:maxres), + & gshieldc_loc_t3(3,-1:maxres), + & gshieldx_t4(3,-1:maxres), gshieldc_t4(3,-1:maxres), + & gshieldc_loc_t4(3,-1:maxres), + & gshieldx_ll(3,-1:maxres), gshieldc_ll(3,-1:maxres), + & gshieldc_loc_ll(3,-1:maxres), & gvdwc_scp(3,maxres),ghpbx(3,maxres),ghpbc(3,maxres), & gloc(maxvar,2),gradcorr(3,maxres),gradxorr(3,maxres), & gradcorr5(3,maxres),gradcorr6(3,maxres), diff --git a/source/cluster/wham/src-M/DIMENSIONS b/source/cluster/wham/src-M/DIMENSIONS index 29e8f24..86f775a 100644 --- a/source/cluster/wham/src-M/DIMENSIONS +++ b/source/cluster/wham/src-M/DIMENSIONS @@ -9,7 +9,7 @@ C Max. number of processors. parameter (maxprocs=16) C Max. number of AA residues integer maxres,maxres2 - parameter (maxres=650) + parameter (maxres=100) C Appr. max. number of interaction sites parameter (maxres2=2*maxres) C Max. number of variables @@ -60,7 +60,7 @@ C Max number of symetric chains parameter (maxperm=120) C Max. number of energy components integer max_ene - parameter (max_ene=24) + parameter (max_ene=25) C Max. number of temperatures integer maxt parameter (maxT=5) diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F index f640679..8203cc4 100644 --- a/source/cluster/wham/src-M/energy_p_new.F +++ b/source/cluster/wham/src-M/energy_p_new.F @@ -22,6 +22,8 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.INTERACT' include 'COMMON.SBRIDGE' include 'COMMON.CHAIN' + include 'COMMON.SHIELD' + include 'COMMON.CONTROL' double precision fact(6) cd write(iout, '(a,i2)')'Calling etotal ipot=',ipot cd print *,'nnt=',nnt,' nct=',nct @@ -47,7 +49,14 @@ C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence). C C Calculate electrostatic (H-bonding) energy of the main chain. C - 106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) + 106 continue + write(iout,*) "shield_mode",shield_mode,ethetacnstr + if (shield_mode.eq.1) then + call set_shield_fac + else if (shield_mode.eq.2) then + call set_shield_fac2 + endif + call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) C C Calculate excluded-volume interaction energy between peptide groups C and side chains. @@ -102,8 +111,20 @@ c print *,ecorr,ecorr5,ecorr6,eturn6 if (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) then call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1) endif -c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t + write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t #ifdef SPLITELE + if (shield_mode.gt.0) then + etot=fact(1)*wsc*(evdw+fact(6)*evdw_t)+fact(1)*wscp*evdw2 + & +welec*fact(1)*ees + & +fact(1)*wvdwpp*evdw1 + & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc + & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 + & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 + & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 + & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d + & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr +C & +wliptran*eliptran + else etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2+welec*fact(1)*ees & +wvdwpp*evdw1 & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc @@ -112,7 +133,20 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr +C & +wliptran*eliptran + endif #else + if (shield_mode.gt.0) then + etot=fact(1)wsc*(evdw+fact(6)*evdw_t)+fact(1)*wscp*evdw2 + & +welec*fact(1)*(ees+evdw1) + & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc + & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 + & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 + & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 + & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d + & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr +C & +wliptran*eliptran + else etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2 & +welec*fact(1)*(ees+evdw1) & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc @@ -121,7 +155,10 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr +C & +wliptran*eliptran + endif #endif + energia(0)=etot energia(1)=evdw #ifdef SCP14 @@ -181,6 +218,7 @@ C #ifdef SPLITELE do i=1,nct do j=1,3 + if (shield_mode.eq.0) then gradc(j,i,icg)=wsc*gvdwc(j,i)+wscp*gvdwc_scp(j,i)+ & welec*fact(1)*gelc(j,i)+wvdwpp*gvdwpp(j,i)+ & wbond*gradb(j,i)+ @@ -193,14 +231,57 @@ C & wcorr6*fact(5)*gradcorr6(j,i)+ & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) + & +wliptran*gliptranc(j,i) gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*fact(2)*gsccorx(j,i) + & +wliptran*gliptranx(j,i) + else + gradc(j,i,icg)=fact(1)*wsc*gvdwc(j,i) + & +fact(1)*wscp*gvdwc_scp(j,i)+ + & welec*fact(1)*gelc(j,i)+fact(1)*wvdwpp*gvdwpp(j,i)+ + & wbond*gradb(j,i)+ + & wstrain*ghpbc(j,i)+ + & wcorr*fact(3)*gradcorr(j,i)+ + & wel_loc*fact(2)*gel_loc(j,i)+ + & wturn3*fact(2)*gcorr3_turn(j,i)+ + & wturn4*fact(3)*gcorr4_turn(j,i)+ + & wcorr5*fact(4)*gradcorr5(j,i)+ + & wcorr6*fact(5)*gradcorr6(j,i)+ + & wturn6*fact(5)*gcorr6_turn(j,i)+ + & wsccor*fact(2)*gsccorc(j,i) + & +wliptran*gliptranc(j,i) + & +welec*gshieldc(j,i) + & +welec*gshieldc_loc(j,i) + & +wcorr*gshieldc_ec(j,i) + & +wcorr*gshieldc_loc_ec(j,i) + & +wturn3*gshieldc_t3(j,i) + & +wturn3*gshieldc_loc_t3(j,i) + & +wturn4*gshieldc_t4(j,i) + & +wturn4*gshieldc_loc_t4(j,i) + & +wel_loc*gshieldc_ll(j,i) + & +wel_loc*gshieldc_loc_ll(j,i) + + gradx(j,i,icg)=fact(1)*wsc*gvdwx(j,i) + & +fact(1)*wscp*gradx_scp(j,i)+ + & wbond*gradbx(j,i)+ + & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ + & wsccor*fact(2)*gsccorx(j,i) + & +wliptran*gliptranx(j,i) + & +welec*gshieldx(j,i) + & +wcorr*gshieldx_ec(j,i) + & +wturn3*gshieldx_t3(j,i) + & +wturn4*gshieldx_t4(j,i) + & +wel_loc*gshieldx_ll(j,i) + + + endif enddo #else - do i=1,nct + do i=1,nct do j=1,3 + if (shield_mode.eq.0) then gradc(j,i,icg)=wsc*gvdwc(j,i)+wscp*gvdwc_scp(j,i)+ & welec*fact(1)*gelc(j,i)+wstrain*ghpbc(j,i)+ & wbond*gradb(j,i)+ @@ -212,11 +293,34 @@ C & wcorr6*fact(5)*gradcorr6(j,i)+ & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) + & +wliptran*gliptranc(j,i) gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*fact(1)*gsccorx(j,i) - enddo + & +wliptran*gliptranx(j,i) + else + gradc(j,i,icg)=fact(1)*wsc*gvdwc(j,i)+ + & fact(1)*wscp*gvdwc_scp(j,i)+ + & welec*fact(1)*gelc(j,i)+wstrain*ghpbc(j,i)+ + & wbond*gradb(j,i)+ + & wcorr*fact(3)*gradcorr(j,i)+ + & wel_loc*fact(2)*gel_loc(j,i)+ + & wturn3*fact(2)*gcorr3_turn(j,i)+ + & wturn4*fact(3)*gcorr4_turn(j,i)+ + & wcorr5*fact(4)*gradcorr5(j,i)+ + & wcorr6*fact(5)*gradcorr6(j,i)+ + & wturn6*fact(5)*gcorr6_turn(j,i)+ + & wsccor*fact(2)*gsccorc(j,i) + & +wliptran*gliptranc(j,i) + gradx(j,i,icg)=fact(1)*wsc*gvdwx(j,i)+ + & fact(1)*wscp*gradx_scp(j,i)+ + & wbond*gradbx(j,i)+ + & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ + & wsccor*fact(1)*gsccorx(j,i) + & +wliptran*gliptranx(j,i) + endif + enddo #endif enddo @@ -774,7 +878,7 @@ C integer icant external icant integer xshift,yshift,zshift - logical energy_dec /.false./ + logical energy_dec /.true./ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 evdw_t=0.0d0 @@ -942,13 +1046,13 @@ c & aux*e2/eps(itypi,itypj) c if (lprn) then sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) epsi=bb(itypi,itypj)**2/aa(itypi,itypj) -c write (iout,'(2(a3,i3,2x),17(0pf7.3))') -c & restyp(itypi),i,restyp(itypj),j, -c & epsi,sigm,chi1,chi2,chip1,chip2, -c & eps1,eps2rt**2,eps3rt**2,sig,sig0ij, -c & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, -c & evdwij -c write (iout,*) "pratial sum", evdw,evdw_t + write (iout,'(2(a3,i3,2x),17(0pf7.3))') + & restyp(itypi),i,restyp(itypj),j, + & epsi,sigm,chi1,chi2,chip1,chip2, + & eps1,eps2rt**2,eps3rt**2,sig,sig0ij, + & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, + & evdwij + write (iout,*) "pratial sum", evdw,evdw_t c endif if (calc_grad) then C Calculate gradient components. @@ -1837,6 +1941,8 @@ C include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' + 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), @@ -1907,15 +2013,15 @@ cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e gcorr_loc(i)=0.0d0 enddo do i=iatel_s,iatel_e - if (i.eq.1) then +C if (i.eq.1) then if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1 - & .or. itype(i+2).eq.ntyp1) cycle - else - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 - & .or. itype(i+2).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 +C & .or. itype(i+2).eq.ntyp1) cycle +C else +C if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 +C & .or. itype(i+2).eq.ntyp1 +C & .or. itype(i-1).eq.ntyp1 &) cycle - endif +C endif if (itel(i).eq.0) goto 1215 dxi=dc(1,i) dyi=dc(2,i) @@ -1935,16 +2041,16 @@ cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e num_conti=0 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) do j=ielstart(i),ielend(i) - if (j.eq.1) then - if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 - & .or.itype(j+2).eq.ntyp1 - &) cycle - else + if (j.le.1) cycle +C if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 +C & .or.itype(j+2).eq.ntyp1 +C &) cycle +C else if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 - & .or.itype(j+2).eq.ntyp1 - & .or.itype(j-1).eq.ntyp1 +C & .or.itype(j+2).eq.ntyp1 +C & .or.itype(j-1).eq.ntyp1 &) cycle - endif +C endif if (itel(j).eq.0) goto 1216 ind=ind+1 iteli=itel(i) @@ -1975,7 +2081,7 @@ C End diagnostics 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 + dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 xj_safe=xj yj_safe=yj zj_safe=zj @@ -1986,7 +2092,7 @@ C End diagnostics 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 + 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 @@ -2032,7 +2138,20 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions c write (iout,*) "i",i,iteli," j",j,itelj," eesij",eesij C 12/26/95 - for the evaluation of multi-body H-bonding interactions ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg) + if (shield_mode.gt.0) then +C fac_shield(i)=0.4 +C fac_shield(j)=0.6 + el1=el1*fac_shield(i)**2*fac_shield(j)**2 + el2=el2*fac_shield(i)**2*fac_shield(j)**2 + eesij=(el1+el2) ees=ees+eesij + else + fac_shield(i)=1.0 + fac_shield(j)=1.0 + eesij=(el1+el2) + ees=ees+eesij + endif +C ees=ees+eesij evdw1=evdw1+evdwij*sss cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, @@ -2055,6 +2174,63 @@ 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 +C enddo +C enddo + 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 + do k=1,3 ghalf=0.5D0*ggg(k) gelc(k,i)=gelc(k,i)+ghalf @@ -2138,15 +2314,19 @@ 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) + & *fac_shield(i)**2*fac_shield(j)**2 enddo do k=1,3 ghalf=0.5D0*ggg(k) gelc(k,i)=gelc(k,i)+ghalf & +(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)+ghalf & +(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 enddo do k=i+1,j-1 do l=1,3 @@ -2394,16 +2574,70 @@ C Contribution to the local-electrostatic energy coming from the i-j pair & +a33*muij(4) cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij cd write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) + 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 C Partial derivatives in virtual-bond dihedral angles gamma if (calc_grad) then + 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 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) + & *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) + & *fac_shield(i)*fac_shield(j) + cd call checkint3(i,j,mu1,mu2,a22,a23,a32,a33,acipa,eel_loc_ij) cd write(iout,*) 'agg ',agg cd write(iout,*) 'aggi ',aggi @@ -2415,6 +2649,8 @@ 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) + & *fac_shield(i)*fac_shield(j) + enddo do k=i+2,j2 do l=1,3 @@ -2425,12 +2661,20 @@ 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) + & *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 ENDIF @@ -2536,9 +2780,21 @@ c fac3=dsqrt(-ael6i)/r0ij**3 fac3=dsqrt(-ael6i)*r3ij ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1) ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2) + 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 c ees0mij=0.0D0 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 @@ -2603,17 +2859,29 @@ C Derivatives due to the contact function 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 C Diagnostics. Comment out or remove after debugging! @@ -2659,6 +2927,8 @@ C Third- and fourth-order contributions from turns include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' + include 'COMMON.SHIELD' + dimension ggg(3) double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2), & e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2), @@ -2667,6 +2937,18 @@ C Third- and fourth-order contributions from turns & aggj(3,4),aggj1(3,4),a_temp(2,2) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2 if (j.eq.i+2) then + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 +C changes suggested by Ana to avoid out of bounds +C & .or.((i+5).gt.nres) +C & .or.((i-1).le.0) +C end of changes suggested by Ana + & .or. itype(i+2).eq.ntyp1 + & .or. itype(i+3).eq.ntyp1 +C & .or. itype(i+5).eq.ntyp1 +C & .or. itype(i).eq.ntyp1 +C & .or. itype(i-1).eq.ntyp1 + & ) goto 179 + CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Third-order contributions @@ -2681,22 +2963,80 @@ cd 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 transpose2(auxmat(1,1),auxmat1(1,1)) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) + 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 eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + eello_t3=0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + cd write (2,*) 'i,',i,' j',j,'eello_turn3', cd & 0.5d0*(pizda(1,1)+pizda(2,2)), cd & ' eello_turn3_num',4*eello_turn3_num if (calc_grad) then +C Derivatives in shield mode + 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)*eello_t3/fac_shield(i) +C & *2.0 + gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,i)*eello_t3/fac_shield(i) + gshieldc_t3(k,iresshield-1)=gshieldc_t3(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)*eello_t3/fac_shield(j) +C & *2.0 + gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,j)*eello_t3/fac_shield(j) + gshieldc_t3(k,iresshield-1)=gshieldc_t3(k,iresshield-1) + & +rlocshield + + enddo + enddo + + do k=1,3 + gshieldc_t3(k,i)=gshieldc_t3(k,i)+ + & grad_shield(k,i)*eello_t3/fac_shield(i) + gshieldc_t3(k,j)=gshieldc_t3(k,j)+ + & grad_shield(k,j)*eello_t3/fac_shield(j) + gshieldc_t3(k,i-1)=gshieldc_t3(k,i-1)+ + & grad_shield(k,i)*eello_t3/fac_shield(i) + gshieldc_t3(k,j-1)=gshieldc_t3(k,j-1)+ + & grad_shield(k,j)*eello_t3/fac_shield(j) + enddo + endif + C Derivatives in gamma(i) call matmat2(EUgder(1,1,i+1),EUg(1,1,i+2),auxmat2(1,1)) call transpose2(auxmat2(1,1),pizda(1,1)) call matmat2(a_temp(1,1),pizda(1,1),pizda(1,1)) gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + C Derivatives in gamma(i+1) call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1)) call transpose2(auxmat2(1,1),pizda(1,1)) call matmat2(a_temp(1,1),pizda(1,1),pizda(1,1)) gel_loc_turn3(i+1)=gel_loc_turn3(i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + C Cartesian derivatives do l=1,3 a_temp(1,1)=aggi(l,1) @@ -2706,6 +3046,8 @@ C Cartesian derivatives call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i)=gcorr3_turn(l,i) & +0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + a_temp(1,1)=aggi1(l,1) a_temp(1,2)=aggi1(l,2) a_temp(2,1)=aggi1(l,3) @@ -2713,6 +3055,8 @@ C Cartesian derivatives call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + a_temp(1,1)=aggj(l,1) a_temp(1,2)=aggj(l,2) a_temp(2,1)=aggj(l,3) @@ -2720,6 +3064,8 @@ C Cartesian derivatives call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j)=gcorr3_turn(l,j) & +0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -2727,9 +3073,24 @@ C Cartesian derivatives call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j1)=gcorr3_turn(l,j1) & +0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + enddo endif + 179 continue else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1) then + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 +C changes suggested by Ana to avoid out of bounds +C & .or.((i+5).gt.nres) +C & .or.((i-1).le.0) +C end of changes suggested by Ana + & .or. itype(i+3).eq.ntyp1 + & .or. itype(i+4).eq.ntyp1 +C & .or. itype(i+5).eq.ntyp1 + & .or. itype(i).eq.ntyp1 +C & .or. itype(i-1).eq.ntyp1 + & ) goto 178 + CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Fourth-order contributions @@ -2757,11 +3118,64 @@ cd call checkint_turn4(i,a_temp,eello_turn4_num) call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) + 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 eello_turn4=eello_turn4-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) + eello_t4=-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) + cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3), cd & ' eello_turn4_num',8*eello_turn4_num C Derivatives in gamma(i) if (calc_grad) then + 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)*eello_t4/fac_shield(i) +C & *2.0 + gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,i)*eello_t4/fac_shield(i) + gshieldc_t4(k,iresshield-1)=gshieldc_t4(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)*eello_t4/fac_shield(j) +C & *2.0 + gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,j)*eello_t4/fac_shield(j) + gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1) + & +rlocshield + + enddo + enddo + + do k=1,3 + gshieldc_t4(k,i)=gshieldc_t4(k,i)+ + & grad_shield(k,i)*eello_t4/fac_shield(i) + gshieldc_t4(k,j)=gshieldc_t4(k,j)+ + & grad_shield(k,j)*eello_t4/fac_shield(j) + gshieldc_t4(k,i-1)=gshieldc_t4(k,i-1)+ + & grad_shield(k,i)*eello_t4/fac_shield(i) + gshieldc_t4(k,j-1)=gshieldc_t4(k,j-1)+ + & grad_shield(k,j)*eello_t4/fac_shield(j) + enddo + endif + 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)) @@ -2769,6 +3183,8 @@ C Derivatives in gamma(i) 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) + & *fac_shield(i)*fac_shield(j) + C Derivatives in gamma(i+1) call transpose2(EUgder(1,1,i+2),e2tder(1,1)) call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1)) @@ -2777,6 +3193,8 @@ C Derivatives in gamma(i+1) call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3) + & *fac_shield(i)*fac_shield(j) + C Derivatives in gamma(i+2) call transpose2(EUgder(1,1,i+3),e3tder(1,1)) call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1)) @@ -2788,6 +3206,8 @@ C Derivatives in gamma(i+2) call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) + C Cartesian derivatives C Derivatives of this turn contributions in DC(i+2) if (j.lt.nres-1) then @@ -2807,6 +3227,8 @@ C Derivatives of this turn contributions in DC(i+2) s3=0.5d0*(pizda(1,1)+pizda(2,2)) ggg(l)=-(s1+s2+s3) gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) + enddo endif C Remaining derivatives of this turn contribution @@ -2825,6 +3247,8 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) + a_temp(1,1)=aggi1(l,1) a_temp(1,2)=aggi1(l,2) a_temp(2,1)=aggi1(l,3) @@ -2839,6 +3263,8 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) + a_temp(1,1)=aggj(l,1) a_temp(1,2)=aggj(l,2) a_temp(2,1)=aggj(l,3) @@ -2853,6 +3279,8 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) + a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -2867,8 +3295,11 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) + enddo endif + 178 continue endif return end @@ -5549,6 +5980,8 @@ c------------------------------------------------------------------------------ include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.SHIELD' + double precision gx(3),gx1(3) logical lprn lprn=.false. @@ -5606,7 +6039,85 @@ C Calculate multi-body contributions to the gradient. & ekont*(coeffp*ees0pij*gacontp_hb3(ll,kk,k)+ & coeffm*ees0mij*gacontm_hb3(ll,kk,k)) enddo - enddo + enddo + if (shield_mode.gt.0) then + j=ees0plist(jj,i) + l=ees0plist(kk,k) +C print *,i,j,fac_shield(i),fac_shield(j), +C &fac_shield(k),fac_shield(l) + if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. + & (fac_shield(k).gt.0).and.(fac_shield(l).gt.0)) then + do ilist=1,ishield_list(i) + iresshield=shield_list(ilist,i) + do m=1,3 + rlocshield=grad_shield_side(m,ilist,i)*ehbcorr/fac_shield(i) +C & *2.0 + gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+ + & rlocshield + & +grad_shield_loc(m,ilist,i)*ehbcorr/fac_shield(i) + gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1) + &+rlocshield + enddo + enddo + do ilist=1,ishield_list(j) + iresshield=shield_list(ilist,j) + do m=1,3 + rlocshield=grad_shield_side(m,ilist,j)*ehbcorr/fac_shield(j) +C & *2.0 + gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+ + & rlocshield + & +grad_shield_loc(m,ilist,j)*ehbcorr/fac_shield(j) + gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1) + & +rlocshield + enddo + enddo + do ilist=1,ishield_list(k) + iresshield=shield_list(ilist,k) + do m=1,3 + rlocshield=grad_shield_side(m,ilist,k)*ehbcorr/fac_shield(k) +C & *2.0 + gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+ + & rlocshield + & +grad_shield_loc(m,ilist,k)*ehbcorr/fac_shield(k) + gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1) + & +rlocshield + enddo + enddo + do ilist=1,ishield_list(l) + iresshield=shield_list(ilist,l) + do m=1,3 + rlocshield=grad_shield_side(m,ilist,l)*ehbcorr/fac_shield(l) +C & *2.0 + gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+ + & rlocshield + & +grad_shield_loc(m,ilist,l)*ehbcorr/fac_shield(l) + gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1) + & +rlocshield + enddo + enddo +C print *,gshieldx(m,iresshield) + do m=1,3 + gshieldc_ec(m,i)=gshieldc_ec(m,i)+ + & grad_shield(m,i)*ehbcorr/fac_shield(i) + gshieldc_ec(m,j)=gshieldc_ec(m,j)+ + & grad_shield(m,j)*ehbcorr/fac_shield(j) + gshieldc_ec(m,i-1)=gshieldc_ec(m,i-1)+ + & grad_shield(m,i)*ehbcorr/fac_shield(i) + gshieldc_ec(m,j-1)=gshieldc_ec(m,j-1)+ + & grad_shield(m,j)*ehbcorr/fac_shield(j) + + gshieldc_ec(m,k)=gshieldc_ec(m,k)+ + & grad_shield(m,k)*ehbcorr/fac_shield(k) + gshieldc_ec(m,l)=gshieldc_ec(m,l)+ + & grad_shield(m,l)*ehbcorr/fac_shield(l) + gshieldc_ec(m,k-1)=gshieldc_ec(m,k-1)+ + & grad_shield(m,k)*ehbcorr/fac_shield(k) + gshieldc_ec(m,l-1)=gshieldc_ec(m,l-1)+ + & grad_shield(m,l)*ehbcorr/fac_shield(l) + + enddo + endif + endif endif ehbcorr=ekont*ees return @@ -7883,4 +8394,340 @@ C----------------------------------------------------------------------- return end C----------------------------------------------------------------------- +C first for shielding is setting of function of side-chains + subroutine set_shield_fac2 + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.IOUNITS' + include 'COMMON.SHIELD' + include 'COMMON.INTERACT' +C this is the squar root 77 devided by 81 the epislion in lipid (in protein) + double precision div77_81/0.974996043d0/, + &div4_81/0.2222222222d0/,sh_frac_dist_grad(3) + +C the vector between center of side_chain and peptide group + double precision pep_side(3),long,side_calf(3), + &pept_group(3),costhet_grad(3),cosphi_grad_long(3), + &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3) +C the line belowe needs to be changed for FGPROC>1 + do i=1,nres-1 + if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle + ishield_list(i)=0 +Cif there two consequtive dummy atoms there is no peptide group between them +C the line below has to be changed for FGPROC>1 + VolumeTotal=0.0 + do k=1,nres + if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle + dist_pep_side=0.0 + dist_side_calf=0.0 + do j=1,3 +C first lets set vector conecting the ithe side-chain with kth side-chain + pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0 +C pep_side(j)=2.0d0 +C and vector conecting the side-chain with its proper calfa + side_calf(j)=c(j,k+nres)-c(j,k) +C side_calf(j)=2.0d0 + pept_group(j)=c(j,i)-c(j,i+1) +C lets have their lenght + dist_pep_side=pep_side(j)**2+dist_pep_side + dist_side_calf=dist_side_calf+side_calf(j)**2 + dist_pept_group=dist_pept_group+pept_group(j)**2 + enddo + dist_pep_side=dsqrt(dist_pep_side) + dist_pept_group=dsqrt(dist_pept_group) + dist_side_calf=dsqrt(dist_side_calf) + do j=1,3 + pep_side_norm(j)=pep_side(j)/dist_pep_side + side_calf_norm(j)=dist_side_calf + enddo +C now sscale fraction + sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield +C print *,buff_shield,"buff" +C now sscale + if (sh_frac_dist.le.0.0) cycle +C If we reach here it means that this side chain reaches the shielding sphere +C Lets add him to the list for gradient + ishield_list(i)=ishield_list(i)+1 +C ishield_list is a list of non 0 side-chain that contribute to factor gradient +C this list is essential otherwise problem would be O3 + shield_list(ishield_list(i),i)=k +C Lets have the sscale value + if (sh_frac_dist.gt.1.0) then + scale_fac_dist=1.0d0 + do j=1,3 + sh_frac_dist_grad(j)=0.0d0 + enddo + else + scale_fac_dist=-sh_frac_dist*sh_frac_dist + & *(2.0d0*sh_frac_dist-3.0d0) + fac_help_scale=6.0d0*(sh_frac_dist-sh_frac_dist**2) + & /dist_pep_side/buff_shield*0.5d0 +C remember for the final gradient multiply sh_frac_dist_grad(j) +C for side_chain by factor -2 ! + do j=1,3 + sh_frac_dist_grad(j)=fac_help_scale*pep_side(j) +C sh_frac_dist_grad(j)=0.0d0 +C scale_fac_dist=1.0d0 +C print *,"jestem",scale_fac_dist,fac_help_scale, +C & sh_frac_dist_grad(j) + enddo + endif +C this is what is now we have the distance scaling now volume... + short=short_r_sidechain(itype(k)) + long=long_r_sidechain(itype(k)) + costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2) + sinthet=short/dist_pep_side*costhet +C now costhet_grad +C costhet=0.6d0 +C sinthet=0.8 + costhet_fac=costhet**3*short**2*(-0.5d0)/dist_pep_side**4 +C sinthet_fac=costhet**2*0.5d0*(short**3/dist_pep_side**4*costhet +C & -short/dist_pep_side**2/costhet) +C costhet_fac=0.0d0 + do j=1,3 + costhet_grad(j)=costhet_fac*pep_side(j) + enddo +C remember for the final gradient multiply costhet_grad(j) +C for side_chain by factor -2 ! +C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1 +C pep_side0pept_group is vector multiplication + pep_side0pept_group=0.0d0 + do j=1,3 + pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j) + enddo + cosalfa=(pep_side0pept_group/ + & (dist_pep_side*dist_side_calf)) + fac_alfa_sin=1.0d0-cosalfa**2 + fac_alfa_sin=dsqrt(fac_alfa_sin) + rkprim=fac_alfa_sin*(long-short)+short +C rkprim=short + +C now costhet_grad + cosphi=1.0d0/dsqrt(1.0d0+rkprim**2/dist_pep_side**2) +C cosphi=0.6 + cosphi_fac=cosphi**3*rkprim**2*(-0.5d0)/dist_pep_side**4 + sinphi=rkprim/dist_pep_side/dsqrt(1.0d0+rkprim**2/ + & dist_pep_side**2) +C sinphi=0.8 + do j=1,3 + cosphi_grad_long(j)=cosphi_fac*pep_side(j) + &+cosphi**3*0.5d0/dist_pep_side**2*(-rkprim) + &*(long-short)/fac_alfa_sin*cosalfa/ + &((dist_pep_side*dist_side_calf))* + &((side_calf(j))-cosalfa* + &((pep_side(j)/dist_pep_side)*dist_side_calf)) +C cosphi_grad_long(j)=0.0d0 + cosphi_grad_loc(j)=cosphi**3*0.5d0/dist_pep_side**2*(-rkprim) + &*(long-short)/fac_alfa_sin*cosalfa + &/((dist_pep_side*dist_side_calf))* + &(pep_side(j)- + &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side) +C cosphi_grad_loc(j)=0.0d0 + enddo +C print *,sinphi,sinthet + VofOverlap=VSolvSphere/2.0d0*(1.0d0-dsqrt(1.0d0-sinphi*sinthet)) + & /VSolvSphere_div +C & *wshield +C now the gradient... + do j=1,3 + grad_shield(j,i)=grad_shield(j,i) +C gradient po skalowaniu + & +(sh_frac_dist_grad(j)*VofOverlap +C gradient po costhet + & +scale_fac_dist*VSolvSphere/VSolvSphere_div/4.0d0* + &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*( + & sinphi/sinthet*costhet*costhet_grad(j) + & +sinthet/sinphi*cosphi*cosphi_grad_long(j))) + & )*wshield +C grad_shield_side is Cbeta sidechain gradient + grad_shield_side(j,ishield_list(i),i)= + & (sh_frac_dist_grad(j)*-2.0d0 + & *VofOverlap + & -scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0* + &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*( + & sinphi/sinthet*costhet*costhet_grad(j) + & +sinthet/sinphi*cosphi*cosphi_grad_long(j))) + & )*wshield + + grad_shield_loc(j,ishield_list(i),i)= + & scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0* + &(1.0d0/(dsqrt(1.0d0-sinphi*sinthet))*( + & sinthet/sinphi*cosphi*cosphi_grad_loc(j) + & )) + & *wshield + enddo + VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist + enddo + fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield) +C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i) + enddo + return + end +C first for shielding is setting of function of side-chains + subroutine set_shield_fac + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.IOUNITS' + include 'COMMON.SHIELD' + include 'COMMON.INTERACT' +C this is the squar root 77 devided by 81 the epislion in lipid (in protein) + double precision div77_81/0.974996043d0/, + &div4_81/0.2222222222d0/,sh_frac_dist_grad(3) + +C the vector between center of side_chain and peptide group + double precision pep_side(3),long,side_calf(3), + &pept_group(3),costhet_grad(3),cosphi_grad_long(3), + &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3) +C the line belowe needs to be changed for FGPROC>1 + do i=1,nres-1 + if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle + ishield_list(i)=0 +Cif there two consequtive dummy atoms there is no peptide group between them +C the line below has to be changed for FGPROC>1 + VolumeTotal=0.0 + do k=1,nres + if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle + dist_pep_side=0.0 + dist_side_calf=0.0 + do j=1,3 +C first lets set vector conecting the ithe side-chain with kth side-chain + pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0 +C pep_side(j)=2.0d0 +C and vector conecting the side-chain with its proper calfa + side_calf(j)=c(j,k+nres)-c(j,k) +C side_calf(j)=2.0d0 + pept_group(j)=c(j,i)-c(j,i+1) +C lets have their lenght + dist_pep_side=pep_side(j)**2+dist_pep_side + dist_side_calf=dist_side_calf+side_calf(j)**2 + dist_pept_group=dist_pept_group+pept_group(j)**2 + enddo + dist_pep_side=dsqrt(dist_pep_side) + dist_pept_group=dsqrt(dist_pept_group) + dist_side_calf=dsqrt(dist_side_calf) + do j=1,3 + pep_side_norm(j)=pep_side(j)/dist_pep_side + side_calf_norm(j)=dist_side_calf + enddo +C now sscale fraction + sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield +C print *,buff_shield,"buff" +C now sscale + if (sh_frac_dist.le.0.0) cycle +C If we reach here it means that this side chain reaches the shielding sphere +C Lets add him to the list for gradient + ishield_list(i)=ishield_list(i)+1 +C ishield_list is a list of non 0 side-chain that contribute to factor gradient +C this list is essential otherwise problem would be O3 + shield_list(ishield_list(i),i)=k +C Lets have the sscale value + if (sh_frac_dist.gt.1.0) then + scale_fac_dist=1.0d0 + do j=1,3 + sh_frac_dist_grad(j)=0.0d0 + enddo + else + scale_fac_dist=-sh_frac_dist*sh_frac_dist + & *(2.0*sh_frac_dist-3.0d0) + fac_help_scale=6.0*(sh_frac_dist-sh_frac_dist**2) + & /dist_pep_side/buff_shield*0.5 +C remember for the final gradient multiply sh_frac_dist_grad(j) +C for side_chain by factor -2 ! + do j=1,3 + sh_frac_dist_grad(j)=fac_help_scale*pep_side(j) +C print *,"jestem",scale_fac_dist,fac_help_scale, +C & sh_frac_dist_grad(j) + enddo + endif +C if ((i.eq.3).and.(k.eq.2)) then +C print *,i,sh_frac_dist,dist_pep,fac_help_scale,scale_fac_dist +C & ,"TU" +C endif + +C this is what is now we have the distance scaling now volume... + short=short_r_sidechain(itype(k)) + long=long_r_sidechain(itype(k)) + costhet=1.0d0/dsqrt(1.0+short**2/dist_pep_side**2) +C now costhet_grad +C costhet=0.0d0 + costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side**4 +C costhet_fac=0.0d0 + do j=1,3 + costhet_grad(j)=costhet_fac*pep_side(j) + enddo +C remember for the final gradient multiply costhet_grad(j) +C for side_chain by factor -2 ! +C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1 +C pep_side0pept_group is vector multiplication + pep_side0pept_group=0.0 + do j=1,3 + pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j) + enddo + cosalfa=(pep_side0pept_group/ + & (dist_pep_side*dist_side_calf)) + fac_alfa_sin=1.0-cosalfa**2 + fac_alfa_sin=dsqrt(fac_alfa_sin) + rkprim=fac_alfa_sin*(long-short)+short +C now costhet_grad + cosphi=1.0d0/dsqrt(1.0+rkprim**2/dist_pep_side**2) + cosphi_fac=cosphi**3*rkprim**2*(-0.5)/dist_pep_side**4 + + do j=1,3 + cosphi_grad_long(j)=cosphi_fac*pep_side(j) + &+cosphi**3*0.5/dist_pep_side**2*(-rkprim) + &*(long-short)/fac_alfa_sin*cosalfa/ + &((dist_pep_side*dist_side_calf))* + &((side_calf(j))-cosalfa* + &((pep_side(j)/dist_pep_side)*dist_side_calf)) + + cosphi_grad_loc(j)=cosphi**3*0.5/dist_pep_side**2*(-rkprim) + &*(long-short)/fac_alfa_sin*cosalfa + &/((dist_pep_side*dist_side_calf))* + &(pep_side(j)- + &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side) + enddo + + VofOverlap=VSolvSphere/2.0d0*(1.0-costhet)*(1.0-cosphi) + & /VSolvSphere_div + & *wshield +C now the gradient... +C grad_shield is gradient of Calfa for peptide groups +C write(iout,*) "shield_compon",i,k,VSolvSphere,scale_fac_dist, +C & costhet,cosphi +C write(iout,*) "cosphi_compon",i,k,pep_side0pept_group, +C & dist_pep_side,dist_side_calf,c(1,k+nres),c(1,k),itype(k) + do j=1,3 + grad_shield(j,i)=grad_shield(j,i) +C gradient po skalowaniu + & +(sh_frac_dist_grad(j) +C gradient po costhet + &-scale_fac_dist*costhet_grad(j)/(1.0-costhet) + &-scale_fac_dist*(cosphi_grad_long(j)) + &/(1.0-cosphi) )*div77_81 + &*VofOverlap +C grad_shield_side is Cbeta sidechain gradient + grad_shield_side(j,ishield_list(i),i)= + & (sh_frac_dist_grad(j)*-2.0d0 + & +scale_fac_dist*costhet_grad(j)*2.0d0/(1.0-costhet) + & +scale_fac_dist*(cosphi_grad_long(j)) + & *2.0d0/(1.0-cosphi)) + & *div77_81*VofOverlap + + grad_shield_loc(j,ishield_list(i),i)= + & scale_fac_dist*cosphi_grad_loc(j) + & *2.0d0/(1.0-cosphi) + & *div77_81*VofOverlap + enddo + VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist + enddo + fac_shield(i)=VolumeTotal*div77_81+div4_81 +C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i) + enddo + return + end +C-------------------------------------------------------------------------- diff --git a/source/cluster/wham/src-M/main_clust.F b/source/cluster/wham/src-M/main_clust.F index 892a6c7..dc77f04 100644 --- a/source/cluster/wham/src-M/main_clust.F +++ b/source/cluster/wham/src-M/main_clust.F @@ -122,7 +122,7 @@ c call flush(iout) ndis=ncon_work*(ncon_work-1)/2 call work_partition(.true.,ndis) #endif - + write(iout,*) "AFTET wort_part",NCON_work DO I=1,NCON_work ICC(I)=I ENDDO @@ -132,6 +132,8 @@ C C CALCULATE DISTANCES C call daread_ccoords(1,ncon_work) + write (iout,*) "AM I HERE" + call flush(iout) ind1=0 DO I=1,NCON_work-1 if (mod(i,100).eq.0) print *,'Calculating RMS i=',i diff --git a/source/cluster/wham/src-M/parmread.F b/source/cluster/wham/src-M/parmread.F index 1368050..4292bd6 100644 --- a/source/cluster/wham/src-M/parmread.F +++ b/source/cluster/wham/src-M/parmread.F @@ -862,9 +862,13 @@ C----------------------- LJK potential -------------------------------- endif goto 50 C---------------------- GB or BP potential ----------------------------- - 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp), - & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp), - & (alp(i),i=1,ntyp) + 30 do i=1,ntyp + read (isidep,*)(eps(i,j),j=i,ntyp) + enddo + read (isidep,*)(sigma0(i),i=1,ntyp) + read (isidep,*)(sigii(i),i=1,ntyp) + read (isidep,*)(chip0(i),i=1,ntyp) + read (isidep,*)(alp(i),i=1,ntyp) C For the GB potential convert sigma'**2 into chi' if (ipot.eq.4) then do i=1,ntyp diff --git a/source/cluster/wham/src-M/probabl.F b/source/cluster/wham/src-M/probabl.F index 0c45402..ec3fbff 100644 --- a/source/cluster/wham/src-M/probabl.F +++ b/source/cluster/wham/src-M/probabl.F @@ -42,7 +42,8 @@ c enddo write (iout,*) me," indstart",indstart(me)," indend",indend(me) call daread_ccoords(indstart(me),indend(me)) #endif -c write (iout,*) "ncon",ncon + write (iout,*) "ncon",ncon + call flush(iout) temper=1.0d0/(beta_h(ib)*1.987D-3) c write (iout,*) "ib",ib," beta_h",beta_h(ib)," temper",temper c quot=1.0d0/(T0*beta_h(ib)*1.987D-3) @@ -54,6 +55,7 @@ c quotl=quotl*quot c kfacl=kfacl*kfac c fT(l)=kfacl/(kfacl-1.0d0+quotl) c enddo +#define DEBUG if (rescale_mode.eq.1) then quot=1.0d0/(T0*beta_h(ib)*1.987D-3) quotl=1.0d0 @@ -88,8 +90,8 @@ c enddo fT(l)=1.12692801104297249644d0/ & dlog(dexp(quotl)+dexp(-quotl)) enddo -c write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft -c call flush(iout) + write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft + call flush(iout) #if defined(FUNCTH) ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ & 320.0d0 @@ -114,25 +116,33 @@ c call flush(iout) do i=1,ncon ii=i #endif -c write (iout,*) "i",i," ii",ii -c call flush(iout) + write (iout,*) "i",i," ii",ii,"ib",ib,scount(me) + call flush(iout) if (ib.eq.1) then do j=1,nres do k=1,3 c(k,j)=allcart(k,j,i) c(k,j+nres)=allcart(k,j+nres,i) + write(iout,*) "coord",i,j,k,allcart(k,j,i),c(k,j), + & c(k,j+nres),allcart(k,j+nres,i) enddo enddo + write(iout,*) "out of j loop" + call flush(iout) do k=1,3 c(k,nres+1)=c(k,1) c(k,nres+nres)=c(k,nres) enddo + write(iout,*) "after nres+nres",nss_all(i) + call flush(iout) nss=nss_all(i) do j=1,nss ihpb(j)=ihpb_all(j,i) jhpb(j)=jhpb_all(j,i) enddo call int_from_cart1(.false.) + write(iout,*) "before etotal" + call flush(iout) call etotal(energia(0),fT) totfree(i)=energia(0) totfree_buf(i)=totfree(i) @@ -207,7 +217,7 @@ c#endif write (iout,*) "evdw2", wscp, evdw2 write (iout,*) "welec", ft(1),welec,ees write (iout,*) "evdw1", wvdwpp,evdw1 - write (iout,*) "ebe" ebe,wang + write (iout,*) "ebe", ebe,wang #endif Fdimless(i)=beta_h(ib)*etot+entfac(ii) Fdimless_buf(i)=Fdimless(i) diff --git a/source/cluster/wham/src-M/read_coords.F b/source/cluster/wham/src-M/read_coords.F index 17f39a3..ac11ac0 100644 --- a/source/cluster/wham/src-M/read_coords.F +++ b/source/cluster/wham/src-M/read_coords.F @@ -178,13 +178,13 @@ c through a ring. #endif endif -#define DEBUG +C#define DEBUG #ifdef DEBUG write (iout,*) "Opening file ",intinname(:ilen(intinname)) write (iout,*) "lenrec",lenrec_in call flush(iout) #endif -#undef DEBUG +C#undef DEBUG c write (iout,*) "maxconf",maxconf i=0 do while (.true.) @@ -283,6 +283,7 @@ c write (iout,*) "nss",nss enddo enddo endif +#define DEBUG #ifdef DEBUG write (iout,'(5hREAD ,i5,3f15.4,i10)') & jj+1,energy(jj+1),entfac(jj+1), @@ -292,6 +293,7 @@ c write (iout,*) "nss",nss write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct) call flush(iout) #endif +#undef DEBUG call add_new_cconf(jjj,jj,jj_old,icount,Next) enddo 101 continue @@ -640,10 +642,11 @@ c------------------------------------------------------------------------------ integer i,j,ij,ii,iii integer len character*16 form,acc - character*32 nam + character*80 nam c c Read conformations off a DA scratchfile. c +#define DEBUG #ifdef DEBUG write (iout,*) "DAREAD_COORDS" write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf @@ -680,7 +683,10 @@ c & nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss), & jhpb_all(i,ij),i=1,nss) call flush(iout) #endif +#undef DEBUG enddo + write (iout,*) "just before leave" + call flush(iout) return end c------------------------------------------------------------------------------ diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F index ea51dbf..6f8f5e4 100644 --- a/source/cluster/wham/src-M/readrtns.F +++ b/source/cluster/wham/src-M/readrtns.F @@ -16,12 +16,13 @@ C include 'COMMON.FREE' include 'COMMON.INTERACT' include "COMMON.SPLITELE" + include 'COMMON.SHIELD' character*320 controlcard,ucase #ifdef MPL include 'COMMON.INFO' #endif integer i,i1,i2,it1,it2 - + double precision pi read (INP,'(a80)') titel call card_concat(controlcard) @@ -36,6 +37,26 @@ C Reading the dimensions of box in x,y,z coordinates c Cutoff range for interactions call reada(controlcard,"R_CUT",r_cut,15.0d0) call reada(controlcard,"LAMBDA",rlamb,0.3d0) +C Shielding mode + call readi(controlcard,'SHIELD',shield_mode,0) + write (iout,*) "SHIELD MODE",shield_mode + if (shield_mode.gt.0) then + pi=3.141592d0 +C VSolvSphere the volume of solving sphere +C print *,pi,"pi" +C rpp(1,1) is the energy r0 for peptide group contact and will be used for it +C there will be no distinction between proline peptide group and normal peptide +C group in case of shielding parameters + VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 + VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 + write (iout,*) VSolvSphere,VSolvSphere_div +C long axis of side chain + do i=1,ntyp + long_r_sidechain(i)=vbldsc0(1,i) + short_r_sidechain(i)=sigma0(i) + enddo + buff_shield=1.0d0 + endif call readi(controlcard,'PDBOUT',outpdb,0) call readi(controlcard,'MOL2OUT',outmol2,0) refstr=(index(controlcard,'REFSTR').gt.0) @@ -103,6 +124,7 @@ C include 'COMMON.CONTACTS' include 'COMMON.TIME1' include 'COMMON.TORCNSTR' + include 'COMMON.SHIELD' #ifdef MPL include 'COMMON.INFO' #endif @@ -118,7 +140,9 @@ C Body C C Read weights of the subsequent energy terms. call card_concat(weightcard) - call reada(weightcard,'WSC',wsc,1.0d0) + write(iout,*) weightcard +C call reada(weightcard,'WSC',wsc,1.0d0) + write(iout,*) wsc call reada(weightcard,'WLONG',wsc,wsc) call reada(weightcard,'WSCP',wscp,1.0d0) call reada(weightcard,'WELEC',welec,1.0D0) @@ -150,6 +174,8 @@ C Read weights of the subsequent energy terms. call reada(weightcard,"V2SS",v2ss,7.61d0) call reada(weightcard,"V3SS",v3ss,13.7d0) call reada(weightcard,"EBR",ebr,-5.50D0) + call reada(weightcard,'WSHIELD',wshield,1.0d0) + write(iout,*) 'WSHIELD',wshield call reada(weightcard,"ATRISS",atriss,0.301D0) call reada(weightcard,"BTRISS",btriss,0.021D0) call reada(weightcard,"CTRISS",ctriss,1.001D0) diff --git a/source/cluster/wham/src-M/work_partition.F b/source/cluster/wham/src-M/work_partition.F index e31db53..ba87bd2 100644 --- a/source/cluster/wham/src-M/work_partition.F +++ b/source/cluster/wham/src-M/work_partition.F @@ -74,12 +74,13 @@ c print *,"N",n," NCON_WORK",ncon_work if (lprint) then write (iout,*) "Partition of work between processors" - do i=0,nprocs-1 - write (iout,'(a,i5,a,i7,a,i7,a,i7)') - & "Processor",i," indstart",indstart(i), - & " indend",indend(i)," count",scount(i) - enddo - endif +C do i=0,nprocs-1 +C write (iout,'(a,i5,a,i7,a,i7,a,i7)') +C & "Processor",i," indstart",indstart(i), +C & " indend",indend(i)," count",scount(i) +C enddo + endif + write(iout,*) "just before leave" return end #endif 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 64288d6..3ad8f52 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -3661,6 +3661,7 @@ C 13-go grudnia roku pamietnego... double precision unmat(3,3) /1.0d0,0.0d0,0.0d0, & 0.0d0,1.0d0,0.0d0, & 0.0d0,0.0d0,1.0d0/ + integer xshift,yshift,zshift c time00=MPI_Wtime() cd write (iout,*) "eelecij",i,j c ind=ind+1 diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index 6c4d182..dcd4d23 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -224,6 +224,11 @@ c detecting NaNQ #ifdef MPL c endif #endif +#define DEBUG +#ifdef DEBUG + call enerprint(energia,fact) +#endif +#undef DEBUG if (calc_grad) then C C Sum up the components of the Cartesian gradient. @@ -265,12 +270,29 @@ C & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) & +wliptran*gliptranc(j,i) + & +welec*gshieldc(j,i) + & +welec*gshieldc_loc(j,i) + & +wcorr*gshieldc_ec(j,i) + & +wcorr*gshieldc_loc_ec(j,i) + & +wturn3*gshieldc_t3(j,i) + & +wturn3*gshieldc_loc_t3(j,i) + & +wturn4*gshieldc_t4(j,i) + & +wturn4*gshieldc_loc_t4(j,i) + & +wel_loc*gshieldc_ll(j,i) + & +wel_loc*gshieldc_loc_ll(j,i) + gradx(j,i,icg)=fact(1)*wsc*gvdwx(j,i) & +fact(1)*wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*fact(2)*gsccorx(j,i) & +wliptran*gliptranx(j,i) + & +welec*gshieldx(j,i) + & +wcorr*gshieldx_ec(j,i) + & +wturn3*gshieldx_t3(j,i) + & +wturn4*gshieldx_t4(j,i) + & +wel_loc*gshieldx_ll(j,i) + endif enddo @@ -309,12 +331,29 @@ C & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) & +wliptran*gliptranc(j,i) + & +welec*gshieldc(j,i) + & +welec*gshieldc_loc(j,i) + & +wcorr*gshieldc_ec(j,i) + & +wcorr*gshieldc_loc_ec(j,i) + & +wturn3*gshieldc_t3(j,i) + & +wturn3*gshieldc_loc_t3(j,i) + & +wturn4*gshieldc_t4(j,i) + & +wturn4*gshieldc_loc_t4(j,i) + & +wel_loc*gshieldc_ll(j,i) + & +wel_loc*gshieldc_loc_ll(j,i) + gradx(j,i,icg)=fact(1)*wsc*gvdwx(j,i)+ & fact(1)*wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*fact(1)*gsccorx(j,i) & +wliptran*gliptranx(j,i) + & +welec*gshieldx(j,i) + & +wcorr*gshieldx_ec(j,i) + & +wturn3*gshieldx_t3(j,i) + & +wturn4*gshieldx_t4(j,i) + & +wel_loc*gshieldx_ll(j,i) + endif enddo #endif @@ -1034,7 +1073,12 @@ C lipbufthick is thickenes of lipid buffore & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 -C write(iout,*) "tu,", i,j,aa_lip(itypi,itypj),bb_lip(itypi,itypj) +C if (aa.ne.aa_aq(itypi,itypj)) then + + write(iout,*) "tu,", i,j,aa_aq(itypi,itypj)-aa, + & bb_aq(itypi,itypj)-bb, + & sslipi,sslipj +C endif C write(iout,*),aa,aa_lip(itypi,itypj),aa_aq(itypi,itypj) C checking the distance @@ -1118,6 +1162,7 @@ c & aux*e2/eps(itypi,itypj) c if (lprn) then sigm=dabs(aa/bb)**(1.0D0/6.0D0) epsi=bb**2/aa +#define DEBUG #ifdef DEBUG write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, @@ -1127,6 +1172,7 @@ c if (lprn) then & evdwij write (iout,*) "partial sum", evdw, evdw_t #endif +#undef DEBUG c endif if (calc_grad) then C Calculate gradient components. @@ -2962,6 +3008,8 @@ C Derivatives due to the contact function & *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) @@ -3121,6 +3169,7 @@ C Derivatives in gamma(i) call transpose2(auxmat2(1,1),pizda(1,1)) call matmat2(a_temp(1,1),pizda(1,1),pizda(1,1)) gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) C Derivatives in gamma(i+1) call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1)) call transpose2(auxmat2(1,1),pizda(1,1)) diff --git a/source/wham/src-M/parmread.F b/source/wham/src-M/parmread.F index 8db78b9..cd2bce9 100644 --- a/source/wham/src-M/parmread.F +++ b/source/wham/src-M/parmread.F @@ -1228,7 +1228,7 @@ c augm(i,j)=0.5D0**(2*expon)*aa(i,j) endif if (lprint) then write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') - & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j), + & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j), & sigma(i,j),r0(i,j),chi(i,j),chi(j,i) endif enddo diff --git a/source/wham/src-M/wham_calc1.F b/source/wham/src-M/wham_calc1.F index eff98da..1d27235 100644 --- a/source/wham/src-M/wham_calc1.F +++ b/source/wham/src-M/wham_calc1.F @@ -1,4 +1,4 @@ - subroutine WHAM_CALC(islice,*) + subroutine WHAM_CALC(islice,*) ! Weighed Histogram Analysis Method (WHAM) code ! Written by A. Liwo based on the work of Kumar et al., ! J.Comput.Chem., 13, 1011 (1992) @@ -36,6 +36,7 @@ c parameter (MaxHdim=200000) include "COMMON.SBRIDGE" include "COMMON.PROT" include "COMMON.ENEPS" + include "COMMON.SHIELD" integer MaxPoint,MaxPointProc parameter (MaxPoint=MaxStr, & MaxPointProc=MaxStr_Proc) @@ -113,6 +114,7 @@ c parameter (MaxHdim=200000) do t=0,MaxN htot(t)=0 enddo +#define DEBUG #ifdef MPI do i=1,scount(me1) #else @@ -820,7 +822,7 @@ C & +ftprim(6)*evdw_t & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ & ftprim(1)*wsccor*esccor ebis=ftbis(1)*wsc*(evdw+ft(6)*evdw_t) - & +ftbis(1)**wscp*evdw2+ + & +ftbis(1)*wscp*evdw2+ & ftbis(1)*welec*ees & +ftbis(1)*wvdwpp*evdw & +ftbis(1)*wtor*etors+ @@ -1284,5 +1286,5 @@ C & +ftprim(6)*evdw_t #endif return - +#undef DEBUG end -- 1.7.9.5