X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-SAXS%2Fenergy_p_new_barrier.F;h=ce803bcbddff8bb5470d26df90b1a0b0d4eaa475;hb=1346fb3959c2eb0a370b11bc6ccad5e4cca27ec9;hp=b61d1e4d6461432c9a4eb0cd8df83934db5ec004;hpb=5836ecdab5a8b95f079bbf6e07374dee3fce8a26;p=unres.git diff --git a/source/unres/src_MD-M-SAXS/energy_p_new_barrier.F b/source/unres/src_MD-M-SAXS/energy_p_new_barrier.F index b61d1e4..ce803bc 100644 --- a/source/unres/src_MD-M-SAXS/energy_p_new_barrier.F +++ b/source/unres/src_MD-M-SAXS/energy_p_new_barrier.F @@ -1,5 +1,5 @@ subroutine etotal(energia) - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifndef ISNAN external proc_proc @@ -10,6 +10,8 @@ cMS$ATTRIBUTES C :: proc_proc #ifdef MPI include "mpif.h" double precision weights_(n_ene) + double precision time00 + integer ierror,ierr #endif include 'COMMON.SETUP' include 'COMMON.IOUNITS' @@ -21,11 +23,19 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.SBRIDGE' include 'COMMON.CHAIN' include 'COMMON.VAR' - include 'COMMON.MD' +c include 'COMMON.MD' + include 'COMMON.QRESTR' include 'COMMON.CONTROL' include 'COMMON.TIME1' include 'COMMON.SPLITELE' include 'COMMON.TORCNSTR' + include 'COMMON.SAXS' + double precision evdw,evdw1,evdw2,evdw2_14,ees,eel_loc, + & eello_turn3,eello_turn4,edfadis,estr,ehpb,ebe,ethetacnstr, + & escloc,etors,edihcnstr,etors_d,esccor,ecorr,ecorr5,ecorr6,eturn6, + & eliptran,Eafmforce,Etube, + & esaxs_constr,ehomology_constr,edfator,edfanei,edfabet + integer n_corr,n_corr1 #ifdef MPI c print*,"ETOTAL Processor",fg_rank," absolute rank",myrank, c & " nfgtasks",nfgtasks @@ -58,6 +68,10 @@ C FG slaves as WEIGHTS array. weights_(21)=wsccor weights_(22)=wtube weights_(26)=wsaxs + weights_(28)=wdfa_dist + weights_(29)=wdfa_tor + weights_(30)=wdfa_nei + weights_(31)=wdfa_beta C FG Master broadcasts the WEIGHTS_ array call MPI_Bcast(weights_(1),n_ene, & MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR) @@ -86,11 +100,21 @@ C FG slaves receive the WEIGHTS array wsccor=weights(21) wtube=weights(22) wsaxs=weights(26) + wdfa_dist=weights_(28) + wdfa_tor=weights_(29) + wdfa_nei=weights_(30) + wdfa_beta=weights_(31) endif time_Bcast=time_Bcast+MPI_Wtime()-time00 time_Bcastw=time_Bcastw+MPI_Wtime()-time00 c call chainbuild_cart endif +#ifndef DFA + edfadis=0.0d0 + edfator=0.0d0 + edfanei=0.0d0 + edfabet=0.0d0 +#endif c print *,'Processor',myrank,' calling etotal ipot=',ipot c print *,'Processor',myrank,' nnt=',nnt,' nct=',nct #else @@ -129,6 +153,32 @@ C C Calculate electrostatic (H-bonding) energy of the main chain. C 107 continue +#ifdef DFA +C BARTEK for dfa test! + if (wdfa_dist.gt.0) then + call edfad(edfadis) + else + edfadis=0 + endif +c print*, 'edfad is finished!', edfadis + if (wdfa_tor.gt.0) then + call edfat(edfator) + else + edfator=0 + endif +c print*, 'edfat is finished!', edfator + if (wdfa_nei.gt.0) then + call edfan(edfanei) + else + edfanei=0 + endif +c print*, 'edfan is finished!', edfanei + if (wdfa_beta.gt.0) then + call edfab(edfabet) + else + edfabet=0 + endif +#endif cmc cmc Sep-06: egb takes care of dynamic ss bonds too cmc @@ -249,6 +299,14 @@ C energy function edihcnstr=0.0d0 if (ndih_constr.gt.0) call etor_constr(edihcnstr) c print *,"Processor",myrank," computed Utor" + if (constr_homology.ge.1) then + call e_modeller(ehomology_constr) +c print *,'iset=',iset,'me=',me,ehomology_constr, +c & 'Processor',fg_rank,' CG group',kolor, +c & ' absolute rank',MyRank + else + ehomology_constr=0.0d0 + endif C C 6/23/01 Calculate double-torsional energy C @@ -266,6 +324,7 @@ C else esccor=0.0d0 endif +#ifdef FOURBODY C print *,"PRZED MULIt" c print *,"Processor",myrank," computed Usccorr" C @@ -294,6 +353,7 @@ c write (iout,*) "MULTIBODY_HB ecorr",ecorr,ecorr5,ecorr6,n_corr, c & n_corr1 c call flush(iout) endif +#endif c print *,"Processor",myrank," computed Ucorr" c write (iout,*) "nsaxs",nsaxs," saxs_mode",saxs_mode if (nsaxs.gt.0 .and. saxs_mode.eq.0) then @@ -389,6 +449,11 @@ C energia(24)=ethetacnstr energia(25)=Etube energia(26)=Esaxs_constr + energia(27)=ehomology_constr + energia(28)=edfadis + energia(29)=edfator + energia(30)=edfanei + energia(31)=edfabet c write (iout,*) "esaxs_constr",energia(26) c Here are the energies showed per procesor if the are more processors c per molecule then we sum it up in sum_energy subroutine @@ -404,7 +469,7 @@ c print *," Processor",myrank," left SUM_ENERGY" end c------------------------------------------------------------------------------- subroutine sum_energy(energia,reduce) - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifndef ISNAN external proc_proc @@ -414,6 +479,8 @@ cMS$ATTRIBUTES C :: proc_proc #endif #ifdef MPI include "mpif.h" + integer ierr + double precision time00 #endif include 'COMMON.SETUP' include 'COMMON.IOUNITS' @@ -427,6 +494,13 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.CONTROL' include 'COMMON.TIME1' logical reduce + integer i + double precision evdw,evdw1,evdw2,evdw2_14,ees,eel_loc, + & eello_turn3,eello_turn4,edfadis,estr,ehpb,ebe,ethetacnstr, + & escloc,etors,edihcnstr,etors_d,esccor,ecorr,ecorr5,ecorr6,eturn6, + & eliptran,Eafmforce,Etube, + & esaxs_constr,ehomology_constr,edfator,edfanei,edfabet + double precision Uconst,etot #ifdef MPI if (nfgtasks.gt.1 .and. reduce) then #ifdef DEBUG @@ -487,6 +561,11 @@ cMS$ATTRIBUTES C :: proc_proc ethetacnstr=energia(24) Etube=energia(25) esaxs_constr=energia(26) + ehomology_constr=energia(27) + edfadis=energia(28) + edfator=energia(29) + edfanei=energia(30) + edfabet=energia(31) #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc @@ -494,7 +573,9 @@ cMS$ATTRIBUTES C :: proc_proc & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & +wbond*estr+wumb*Uconst+wsccor*esccor+wliptran*eliptran+Eafmforce - & +ethetacnstr+wtube*Etube+wsaxs*esaxs_constr + & +ethetacnstr+wtube*Etube+wsaxs*esaxs_constr+ehomology_constr + & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei + & +wdfa_beta*edfabet #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc @@ -503,7 +584,9 @@ cMS$ATTRIBUTES C :: proc_proc & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & +wbond*estr+wumb*Uconst+wsccor*esccor+wliptran*eliptran & +Eafmforce - & +ethetacnstr+wtube*Etube+wsaxs*esaxs_constr + & +ethetacnstr+wtube*Etube+wsaxs*esaxs_constr+ehomology_constr + & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei + & +wdfa_beta*edfabet #endif energia(0)=etot c detecting NaNQ @@ -529,7 +612,7 @@ c detecting NaNQ end c------------------------------------------------------------------------------- subroutine sum_gradient - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifndef ISNAN external proc_proc @@ -539,6 +622,8 @@ cMS$ATTRIBUTES C :: proc_proc #endif #ifdef MPI include 'mpif.h' + integer ierror,ierr + double precision time00,time01 #endif double precision gradbufc(3,-1:maxres),gradbufx(3,-1:maxres), & glocbuf(4*maxres),gradbufc_sum(3,-1:maxres) @@ -555,6 +640,16 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.TIME1' include 'COMMON.MAXGRAD' include 'COMMON.SCCOR' +c include 'COMMON.MD' + include 'COMMON.QRESTR' + integer i,j,k + double precision scalar + double precision gvdwc_norm,gvdwc_scp_norm,gelc_norm,gvdwpp_norm, + &gradb_norm,ghpbc_norm,gradcorr_norm,gel_loc_norm,gcorr3_turn_norm, + &gcorr4_turn_norm,gradcorr5_norm,gradcorr6_norm, + &gcorr6_turn_norm,gsccorrc_norm,gscloc_norm,gvdwx_norm, + &gradx_scp_norm,ghpbx_norm,gradxorr_norm,gsccorrx_norm, + &gsclocx_norm #ifdef TIMING time01=MPI_Wtime() #endif @@ -623,16 +718,13 @@ c enddo & wstrain*ghpbc(j,i) & +wliptran*gliptranc(j,i) & +gradafm(j,i) - & +welec*gshieldc(j,i) - & +wcorr*gshieldc_ec(j,i) - & +wturn3*gshieldc_t3(j,i) - & +wturn4*gshieldc_t4(j,i) - & +wel_loc*gshieldc_ll(j,i) + & +welec*gshieldc(j,i) + & +wcorr*gshieldc_ec(j,i) + & +wturn3*gshieldc_t3(j,i) + & +wturn4*gshieldc_t4(j,i) + & +wel_loc*gshieldc_ll(j,i) & +wtube*gg_tube(j,i) & +wsaxs*gsaxsc(j,i) - - - enddo enddo #else @@ -656,12 +748,18 @@ c enddo & +wel_loc*gshieldc_ll(j,i) & +wtube*gg_tube(j,i) & +wsaxs*gsaxsc(j,i) - - - enddo enddo #endif + do i=1,nct + do j=1,3 + gradbufc(j,i)=gradbufc(j,i)+ + & wdfa_dist*gdfad(j,i)+ + & wdfa_tor*gdfat(j,i)+ + & wdfa_nei*gdfan(j,i)+ + & wdfa_beta*gdfab(j,i) + enddo + enddo #ifdef DEBUG write (iout,*) "gradc from gradbufc" do i=1,nres @@ -875,6 +973,14 @@ C print *,gradafm(1,13),"AFM" enddo enddo + if (constr_homology.gt.0) then + do i=1,nct + do j=1,3 + gradc(j,i,icg)=gradc(j,i,icg)+duscdiff(j,i) + gradx(j,i,icg)=gradx(j,i,icg)+duscdiffx(j,i) + enddo + enddo + endif #ifdef DEBUG write (iout,*) "gradc gradx gloc after adding" do i=1,nres @@ -984,13 +1090,13 @@ c gradcorr5_max=0.0d0 gradcorr6_max=0.0d0 gcorr6_turn_max=0.0d0 - gsccorc_max=0.0d0 + gsccorrc_max=0.0d0 gscloc_max=0.0d0 gvdwx_max=0.0d0 gradx_scp_max=0.0d0 ghpbx_max=0.0d0 gradxorr_max=0.0d0 - gsccorx_max=0.0d0 + gsccorrx_max=0.0d0 gsclocx_max=0.0d0 do i=1,nct gvdwc_norm=dsqrt(scalar(gvdwc(1,i),gvdwc(1,i))) @@ -1022,13 +1128,13 @@ c if (gradcorr5_norm.gt.gradcorr5_max) & gradcorr5_max=gradcorr5_norm gradcorr6_norm=dsqrt(scalar(gradcorr6(1,i),gradcorr6(1,i))) - if (gradcorr6_norm.gt.gradcorr6_max) gcorr6_max=gradcorr6_norm + if (gradcorr6_norm.gt.gradcorr6_max)gradcorr6_max=gradcorr6_norm gcorr6_turn_norm=dsqrt(scalar(gcorr6_turn(1,i), & gcorr6_turn(1,i))) if (gcorr6_turn_norm.gt.gcorr6_turn_max) & gcorr6_turn_max=gcorr6_turn_norm - gsccorr_norm=dsqrt(scalar(gsccorc(1,i),gsccorc(1,i))) - if (gsccorr_norm.gt.gsccorr_max) gsccorr_max=gsccorr_norm + gsccorrc_norm=dsqrt(scalar(gsccorc(1,i),gsccorc(1,i))) + if (gsccorrc_norm.gt.gsccorrc_max) gsccorrc_max=gsccorrc_norm gscloc_norm=dsqrt(scalar(gscloc(1,i),gscloc(1,i))) if (gscloc_norm.gt.gscloc_max) gscloc_max=gscloc_norm gvdwx_norm=dsqrt(scalar(gvdwx(1,i),gvdwx(1,i))) @@ -1054,9 +1160,9 @@ c write (istat,'(1h#,21f10.2)') gvdwc_max,gvdwc_scp_max, & gelc_max,gvdwpp_max,gradb_max,ghpbc_max, & gradcorr_max,gel_loc_max,gcorr3_turn_max,gcorr4_turn_max, - & gradcorr5_max,gradcorr6_max,gcorr6_turn_max,gsccorc_max, + & gradcorr5_max,gradcorr6_max,gcorr6_turn_max,gsccorrc_max, & gscloc_max,gvdwx_max,gradx_scp_max,ghpbx_max,gradxorr_max, - & gsccorx_max,gsclocx_max + & gsccorrx_max,gsclocx_max close(istat) if (gvdwc_max.gt.1.0d4) then write (iout,*) "gvdwc gvdwx gradb gradbx" @@ -1083,12 +1189,18 @@ c end c------------------------------------------------------------------------------- subroutine rescale_weights(t_bath) - implicit real*8 (a-h,o-z) + implicit none +#ifdef MPI + include 'mpif.h' + integer ierror +#endif include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.FFIELD' include 'COMMON.SBRIDGE' include 'COMMON.CONTROL' + double precision t_bath + double precision facT,facT2,facT3,facT4,facT5 double precision kfac /2.4d0/ double precision x,x2,x3,x4,x5,licznik /1.12692801104297249644/ c facT=temp0/t_bath @@ -1148,13 +1260,19 @@ c write (iout,*) "t_bath",t_bath," temp0",temp0," wumb",wumb end C------------------------------------------------------------------------ subroutine enerprint(energia) - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.FFIELD' include 'COMMON.SBRIDGE' - include 'COMMON.MD' + include 'COMMON.QRESTR' double precision energia(0:n_ene) + double precision evdw,evdw1,evdw2,evdw2_14,ees,eel_loc, + & eello_turn3,eello_turn4,edfadis,estr,ehpb,ebe,ethetacnstr, + & escloc,etors,edihcnstr,etors_d,esccor,ecorr,ecorr5,ecorr6, + & eello_turn6, + & eliptran,Eafmforce,Etube, + & esaxs,ehomology_constr,edfator,edfanei,edfabet,etot etot=energia(0) evdw=energia(1) evdw2=energia(2) @@ -1188,83 +1306,125 @@ C------------------------------------------------------------------------ ethetacnstr=energia(24) etube=energia(25) esaxs=energia(26) + ehomology_constr=energia(27) +C Bartek + edfadis = energia(28) + edfator = energia(29) + edfanei = energia(30) + edfabet = energia(31) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp, & estr,wbond,ebe,wang, & escloc,wscloc,etors,wtor,etors_d,wtor_d,ehpb,wstrain, +#ifdef FOURBODY & ecorr,wcorr, - & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, - & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr, - & ethetacnstr,ebr*nss,Uconst,wumb,eliptran,wliptran,Eafmforc, - & etube,wtube,esaxs,wsaxs, + & ecorr5,wcorr5,ecorr6,wcorr6, +#endif + & eel_loc,wel_loc,eello_turn3,wturn3, + & eello_turn4,wturn4, +#ifdef FOURBODY + & eello_turn6,wturn6, +#endif + & esccor,wsccor,edihcnstr, + & ethetacnstr,ebr*nss,Uconst,wumb,eliptran,wliptran,Eafmforce, + & etube,wtube,esaxs,wsaxs,ehomology_constr, + & edfadis,wdfa_dist,edfator,wdfa_tor,edfanei,wdfa_nei, + & edfabet,wdfa_beta, & etot 10 format (/'Virtual-chain energies:'// - & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ - & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ - & 'EES= ',1pE16.6,' WEIGHT=',1pD16.6,' (p-p)'/ - & 'EVDWPP=',1pE16.6,' WEIGHT=',1pD16.6,' (p-p VDW)'/ - & 'ESTR= ',1pE16.6,' WEIGHT=',1pD16.6,' (stretching)'/ - & 'EBE= ',1pE16.6,' WEIGHT=',1pD16.6,' (bending)'/ - & 'ESC= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC local)'/ - & 'ETORS= ',1pE16.6,' WEIGHT=',1pD16.6,' (torsional)'/ - & 'ETORSD=',1pE16.6,' WEIGHT=',1pD16.6,' (double torsional)'/ - & 'EHBP= ',1pE16.6,' WEIGHT=',1pD16.6, + & 'EVDW= ',1pE16.6,' WEIGHT=',1pE16.6,' (SC-SC)'/ + & 'EVDW2= ',1pE16.6,' WEIGHT=',1pE16.6,' (SC-p)'/ + & 'EES= ',1pE16.6,' WEIGHT=',1pE16.6,' (p-p)'/ + & 'EVDWPP=',1pE16.6,' WEIGHT=',1pE16.6,' (p-p VDW)'/ + & 'ESTR= ',1pE16.6,' WEIGHT=',1pE16.6,' (stretching)'/ + & 'EBE= ',1pE16.6,' WEIGHT=',1pE16.6,' (bending)'/ + & 'ESC= ',1pE16.6,' WEIGHT=',1pE16.6,' (SC local)'/ + & 'ETORS= ',1pE16.6,' WEIGHT=',1pE16.6,' (torsional)'/ + & 'ETORSD=',1pE16.6,' WEIGHT=',1pE16.6,' (double torsional)'/ + & 'EHBP= ',1pE16.6,' WEIGHT=',1pE16.6, & ' (SS bridges & dist. cnstr.)'/ - & 'ECORR4=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ - & 'ECORR5=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ - & 'ECORR6=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ - & 'EELLO= ',1pE16.6,' WEIGHT=',1pD16.6,' (electrostatic-local)'/ - & 'ETURN3=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 3rd order)'/ - & 'ETURN4=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 4th order)'/ - & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ - & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ +#ifdef FOURBODY + & 'ECORR4=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ + & 'ECORR5=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ + & 'ECORR6=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ +#endif + & 'EELLO= ',1pE16.6,' WEIGHT=',1pE16.6,' (electrostatic-local)'/ + & 'ETURN3=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 3rd order)'/ + & 'ETURN4=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 4th order)'/ +#ifdef FOURBODY + & 'ETURN6=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 6th order)'/ +#endif + & 'ESCCOR=',1pE16.6,' WEIGHT=',1pE16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (virtual-bond dihedral angle restraints)'/ & 'ETHETC=',1pE16.6,' (virtual-bond angle restraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ - & 'UCONST=',1pE16.6,' WEIGHT=',1pD16.6' (umbrella restraints)'/ - & 'ELT= ',1pE16.6,' WEIGHT=',1pD16.6,' (Lipid transfer)'/ + & 'UCONST=',1pE16.6,' WEIGHT=',1pE16.6' (umbrella restraints)'/ + & 'ELT= ',1pE16.6,' WEIGHT=',1pE16.6,' (Lipid transfer)'/ & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ - & 'ETUBE= ',1pE16.6,' WEIGHT=',1pD16.6,' (tube confinment)'/ - & 'E_SAXS=',1pE16.6,' WEIGHT=',1pD16.6,' (SAXS restraints)'/ + & 'ETUBE= ',1pE16.6,' WEIGHT=',1pE16.6,' (tube confinment)'/ + & 'E_SAXS=',1pE16.6,' WEIGHT=',1pE16.6,' (SAXS restraints)'/ + & 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/ + & 'EDFAD= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA distance energy)'/ + & 'EDFAT= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA torsion energy)'/ + & 'EDFAN= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA NCa energy)'/ + & 'EDFAB= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA Beta energy)'/ & 'ETOT= ',1pE16.6,' (total)') #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec, & estr,wbond,ebe,wang, & escloc,wscloc,etors,wtor,etors_d,wtor_d,ehpb,wstrain, +#ifdef FOURBODY & ecorr,wcorr, - & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, - & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr, + & ecorr5,wcorr5,ecorr6,wcorr6, +#endif + & eel_loc,wel_loc,eello_turn3,wturn3, + & eello_turn4,wturn4, +#ifdef FOURBODY + & eello_turn6,wturn6, +#endif + & esccor,wsccor,edihcnstr, & ethetacnstr,ebr*nss,Uconst,wumb,eliptran,wliptran,Eafmforc, - & etube,wtube,esaxs,wsaxs, + & etube,wtube,esaxs,wsaxs,ehomology_constr, + & edfadis,wdfa_dist,edfator,wdfa_tor,edfanei,wdfa_nei, + & edfabet,wdfa_beta, & etot 10 format (/'Virtual-chain energies:'// - & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ - & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ - & 'EES= ',1pE16.6,' WEIGHT=',1pD16.6,' (p-p)'/ - & 'ESTR= ',1pE16.6,' WEIGHT=',1pD16.6,' (stretching)'/ - & 'EBE= ',1pE16.6,' WEIGHT=',1pD16.6,' (bending)'/ - & 'ESC= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC local)'/ - & 'ETORS= ',1pE16.6,' WEIGHT=',1pD16.6,' (torsional)'/ - & 'ETORSD=',1pE16.6,' WEIGHT=',1pD16.6,' (double torsional)'/ - & 'EHBP= ',1pE16.6,' WEIGHT=',1pD16.6, + & 'EVDW= ',1pE16.6,' WEIGHT=',1pE16.6,' (SC-SC)'/ + & 'EVDW2= ',1pE16.6,' WEIGHT=',1pE16.6,' (SC-p)'/ + & 'EES= ',1pE16.6,' WEIGHT=',1pE16.6,' (p-p)'/ + & 'ESTR= ',1pE16.6,' WEIGHT=',1pE16.6,' (stretching)'/ + & 'EBE= ',1pE16.6,' WEIGHT=',1pE16.6,' (bending)'/ + & 'ESC= ',1pE16.6,' WEIGHT=',1pE16.6,' (SC local)'/ + & 'ETORS= ',1pE16.6,' WEIGHT=',1pE16.6,' (torsional)'/ + & 'ETORSD=',1pE16.6,' WEIGHT=',1pE16.6,' (double torsional)'/ + & 'EHBP= ',1pE16.6,' WEIGHT=',1pE16.6, & ' (SS bridges & dist. restr.)'/ - & 'ECORR4=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ - & 'ECORR5=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ - & 'ECORR6=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ - & 'EELLO= ',1pE16.6,' WEIGHT=',1pD16.6,' (electrostatic-local)'/ - & 'ETURN3=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 3rd order)'/ - & 'ETURN4=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 4th order)'/ - & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ - & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ +#ifdef FOURBODY + & 'ECORR4=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ + & 'ECORR5=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ + & 'ECORR6=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ +#endif + & 'EELLO= ',1pE16.6,' WEIGHT=',1pE16.6,' (electrostatic-local)'/ + & 'ETURN3=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 3rd order)'/ + & 'ETURN4=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 4th order)'/ +#ifdef FOURBODY + & 'ETURN6=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 6th order)'/ +#endif + & 'ESCCOR=',1pE16.6,' WEIGHT=',1pE16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (virtual-bond dihedral angle restraints)'/ & 'ETHETC=',1pE16.6,' (virtual-bond angle restraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ - & 'UCONST=',1pE16.6,' WEIGHT=',1pD16.6' (umbrella restraints)'/ - & 'ELT= ',1pE16.6,' WEIGHT=',1pD16.6,' (Lipid transfer)'/ + & 'UCONST=',1pE16.6,' WEIGHT=',1pE16.6' (umbrella restraints)'/ + & 'ELT= ',1pE16.6,' WEIGHT=',1pE16.6,' (Lipid transfer)'/ & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ - & 'ETUBE= ',1pE16.6,' WEIGHT=',1pD16.6,' (tube confinment)'/ - & 'E_SAXS=',1pE16.6,' WEIGHT=',1pD16.6,' (SAXS restraints)'/ + & 'ETUBE= ',1pE16.6,' WEIGHT=',1pE16.6,' (tube confinment)'/ + & 'E_SAXS=',1pE16.6,' WEIGHT=',1pE16.6,' (SAXS restraints)'/ + & 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/ + & 'EDFAD= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA distance energy)'/ + & 'EDFAT= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA torsion energy)'/ + & 'EDFAN= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA NCa energy)'/ + & 'EDFAB= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA Beta energy)'/ & 'ETOT= ',1pE16.6,' (total)') #endif return @@ -1275,7 +1435,8 @@ C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the LJ potential of interaction. C - implicit real*8 (a-h,o-z) + implicit none + double precision accur include 'DIMENSIONS' parameter (accur=1.0d-10) include 'COMMON.GEO' @@ -1288,8 +1449,18 @@ C include 'COMMON.SBRIDGE' include 'COMMON.NAMES' include 'COMMON.IOUNITS' + include 'COMMON.SPLITELE' +#ifdef FOURBODY include 'COMMON.CONTACTS' - dimension gg(3) + include 'COMMON.CONTMAT' +#endif + double precision gg(3) + double precision evdw,evdwij + integer i,j,k,itypi,itypj,itypi1,num_conti,iint + double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, + & sigij,r0ij,rcut,sqrij,sss1,sssgrad1 + double precision fcont,fprimcont + double precision sscale,sscagrad c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e @@ -1316,6 +1487,11 @@ cd & 'iend=',iend(i,iint) C Change 12/1/95 to calculate four-body interactions rij=xj*xj+yj*yj+zj*zj rrij=1.0D0/rij + sqrij=dsqrt(rij) + sss1=sscale(sqrij,r_cut_int) + if (sss1.eq.0.0d0) cycle + sssgrad1=sscagrad(sqrij,r_cut_int) + c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj eps0ij=eps(itypi,itypj) fac=rrij**expon2 @@ -1329,11 +1505,12 @@ cd write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)') cd & restyp(itypi),i,restyp(itypj),j,a(itypi,itypj), cd & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm, cd & (c(k,i),k=1,3),(c(k,j),k=1,3) - evdw=evdw+evdwij + evdw=evdw+sss1*evdwij C C Calculate the components of the gradient in DC and X C fac=-rrij*(e1+evdwij) + & +evdwij*sssgrad1/sqrij gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac @@ -1349,6 +1526,7 @@ cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l) cgrad enddo cgrad enddo C +#ifdef FOURBODY C 12/1/95, revised on 5/20/97 C C Calculate the contact function. The ith column of the array JCONT will @@ -1404,10 +1582,13 @@ cd write (iout,'(2i3,3f10.5)') cd & i,j,(gacont(kk,num_conti,i),kk=1,3) endif endif +#endif enddo ! j enddo ! iint C Change 12/1/95 +#ifdef FOURBODY num_cont(i)=num_conti +#endif enddo ! i do i=1,nct do j=1,3 @@ -1432,7 +1613,7 @@ C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the LJK potential of interaction. C - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -1442,8 +1623,14 @@ C include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' - dimension gg(3) + include 'COMMON.SPLITELE' + double precision gg(3) + double precision evdw,evdwij + integer i,j,k,itypi,itypj,itypi1,iint + double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, + & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1 logical scheck + double precision sscale,sscagrad c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e @@ -1468,6 +1655,9 @@ C e_augm=augm(itypi,itypj)*fac_augm r_inv_ij=dsqrt(rrij) rij=1.0D0/r_inv_ij + sss1=sscale(rij,r_cut_int) + if (sss1.eq.0.0d0) cycle + sssgrad1=sscagrad(rij,r_cut_int) r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon C have you changed here? @@ -1486,6 +1676,7 @@ C C Calculate the components of the gradient in DC and X C fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2) + & +evdwij*sssgrad1*r_inv_ij gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac @@ -1517,7 +1708,7 @@ C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the Berne-Pechukas potential of interaction. C - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -1528,7 +1719,14 @@ C include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.CALC' + include 'COMMON.SPLITELE' + integer icall common /srutu/ icall + double precision evdw + integer itypi,itypj,itypi1,iint,ind + double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi, + & sss1,sssgrad1 + double precision sscale,sscagrad c double precision rrsave(maxdim) logical lprn evdw=0.0D0 @@ -1594,6 +1792,9 @@ cd else cd rrij=rrsave(ind) cd endif rij=dsqrt(rrij) + sss1=sscale(1.0d0/rij,r_cut_int) + if (sss1.eq.0.0d0) cycle + sssgrad1=sscagrad(1.0d0/rij,r_cut_int) C Calculate the angle-dependent terms of energy & contributions to derivatives. call sc_angular C Calculate whole angle-dependent part of epsilon and contributions @@ -1622,13 +1823,14 @@ C Calculate gradient components. fac=-expon*(e1+evdwij) sigder=fac/sigsq fac=rrij*fac + & +evdwij*sssgrad1*rij C Calculate radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac C Calculate the angular part of the gradient and sum add the contributions C to the appropriate components of the Cartesian gradient. - call sc_grad + call sc_grad_scale(sss1) enddo ! j enddo ! iint enddo ! i @@ -1641,7 +1843,7 @@ C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the Gay-Berne potential of interaction. C - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -1656,8 +1858,14 @@ C include 'COMMON.SPLITELE' include 'COMMON.SBRIDGE' logical lprn - integer xshift,yshift,zshift - + integer xshift,yshift,zshift,subchap + double precision evdw + integer itypi,itypj,itypi1,iint,ind + double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi + double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij, + & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe, + & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip + double precision dist,sscale,sscagrad,sscagradlip,sscalelip evdw=0.0D0 ccccc energy_dec=.false. C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -1920,12 +2128,11 @@ c write (iout,*) "j",j," dc_norm", c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) - sss=sscale((1.0d0/rij)/sigma(itypi,itypj)) - sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj)) - + sss=sscale(1.0d0/rij,r_cut_int) c write (iout,'(a7,4f8.3)') c & "ssscale",sss,((1.0d0/rij)/sigma(itypi,itypj)),r_cut,rlamb - if (sss.gt.0.0d0) then + if (sss.eq.0.0d0) cycle + sssgrad=sscagrad(1.0d0/rij,r_cut_int) C Calculate angle-dependent terms of energy and contributions to their C derivatives. call sc_angular @@ -1972,8 +2179,8 @@ c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 & evdwij endif - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') - & 'evdw',i,j,evdwij + if (energy_dec) write (iout,'(a,2i5,3f10.5)') + & 'r sss evdw',i,j,rij,sss,evdwij C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 @@ -1982,13 +2189,13 @@ C Calculate gradient components. fac=rij*fac c print '(2i4,6f8.4)',i,j,sss,sssgrad* c & evdwij,fac,sigma(itypi,itypj),expon - fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij + fac=fac+evdwij*sssgrad*rij c fac=0.0d0 C Calculate the radial part of the gradient gg_lipi(3)=eps1*(eps2rt*eps2rt) - &*(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* - & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) - &+faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) gg_lipj(3)=ssgradlipj*gg_lipi(3) gg_lipi(3)=gg_lipi(3)*ssgradlipi C gg_lipi(3)=0.0d0 @@ -1997,8 +2204,7 @@ C gg_lipj(3)=0.0d0 gg(2)=yj*fac gg(3)=zj*fac C Calculate angular part of the gradient. - call sc_grad - endif + call sc_grad_scale(sss) ENDIF ! dyn_ss enddo ! j enddo ! iint @@ -2016,7 +2222,7 @@ C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the Gay-Berne-Vorobjev potential of interaction. C - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -2027,9 +2233,19 @@ C include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.CALC' - integer xshift,yshift,zshift + include 'COMMON.SPLITELE' + integer xshift,yshift,zshift,subchap + integer icall common /srutu/ icall logical lprn + double precision evdw + integer itypi,itypj,itypi1,iint,ind + double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij, + & xi,yi,zi,fac_augm,e_augm + double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij, + & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe, + & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip,sssgrad1 + double precision dist,sscale,sscagrad,sscagradlip,sscalelip evdw=0.0D0 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 @@ -2188,6 +2404,9 @@ C write(iout,*) "tu,", i,j,aa,bb,aa_lip(itypi,itypj),sslipi,sslipj dzj=dc_norm(3,nres+j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) + sss=sscale(1.0d0/rij,r_cut_int) + if (sss.eq.0.0d0) cycle + sssgrad=sscagrad(1.0d0/rij,r_cut_int) C Calculate angle-dependent terms of energy and contributions to their C derivatives. call sc_angular @@ -2228,13 +2447,13 @@ C Calculate gradient components. fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder fac=rij*fac-2*expon*rrij*e_augm - fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij + fac=fac+(evdwij+e_augm)*sssgrad*rij C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac C Calculate angular part of the gradient. - call sc_grad + call sc_grad_scale(sss) enddo ! j enddo ! iint enddo ! i @@ -2383,7 +2602,7 @@ C include 'COMMON.SBRIDGE' include 'COMMON.NAMES' include 'COMMON.IOUNITS' - include 'COMMON.CONTACTS' +c include 'COMMON.CONTACTS' dimension gg(3) cd print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct evdw=0.0D0 @@ -2457,7 +2676,7 @@ C include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' - include 'COMMON.CONTACTS' +c include 'COMMON.CONTACTS' include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' @@ -2538,8 +2757,8 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) zj=zj_safe-zmedi endif rij=xj*xj+yj*yj+zj*zj - sss=sscale(sqrt(rij)) - sssgrad=sscagrad(sqrt(rij)) + sss=sscale(sqrt(rij),r_cut_int) + sssgrad=sscagrad(sqrt(rij),r_cut_int) if (rij.lt.r0ijsq) then evdw1ij=0.25d0*(rij-r0ijsq)**2 fac=rij-r0ijsq @@ -2764,90 +2983,6 @@ c & " ivec_count",(ivec_count(i),i=0,nfgtasks1-1) #endif return end -C----------------------------------------------------------------------------- - subroutine check_vecgrad - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' - include 'COMMON.IOUNITS' - include 'COMMON.GEO' - include 'COMMON.VAR' - include 'COMMON.LOCAL' - include 'COMMON.CHAIN' - include 'COMMON.VECTORS' - dimension uygradt(3,3,2,maxres),uzgradt(3,3,2,maxres) - dimension uyt(3,maxres),uzt(3,maxres) - dimension uygradn(3,3,2),uzgradn(3,3,2),erij(3) - double precision delta /1.0d-7/ - call vec_and_deriv -cd do i=1,nres -crc write(iout,'(2i5,2(3f10.5,5x))') i,1,dc_norm(:,i) -crc write(iout,'(2i5,2(3f10.5,5x))') i,2,uy(:,i) -crc write(iout,'(2i5,2(3f10.5,5x)/)')i,3,uz(:,i) -cd write(iout,'(2i5,2(3f10.5,5x))') i,1, -cd & (dc_norm(if90,i),if90=1,3) -cd write(iout,'(2i5,2(3f10.5,5x))') i,2,(uy(if90,i),if90=1,3) -cd write(iout,'(2i5,2(3f10.5,5x)/)')i,3,(uz(if90,i),if90=1,3) -cd write(iout,'(a)') -cd enddo - do i=1,nres - do j=1,2 - do k=1,3 - do l=1,3 - uygradt(l,k,j,i)=uygrad(l,k,j,i) - uzgradt(l,k,j,i)=uzgrad(l,k,j,i) - enddo - enddo - enddo - enddo - call vec_and_deriv - do i=1,nres - do j=1,3 - uyt(j,i)=uy(j,i) - uzt(j,i)=uz(j,i) - enddo - enddo - do i=1,nres -cd write (iout,*) 'i=',i - do k=1,3 - erij(k)=dc_norm(k,i) - enddo - do j=1,3 - do k=1,3 - dc_norm(k,i)=erij(k) - enddo - dc_norm(j,i)=dc_norm(j,i)+delta -c fac=dsqrt(scalar(dc_norm(1,i),dc_norm(1,i))) -c do k=1,3 -c dc_norm(k,i)=dc_norm(k,i)/fac -c enddo -c write (iout,*) (dc_norm(k,i),k=1,3) -c write (iout,*) (erij(k),k=1,3) - call vec_and_deriv - do k=1,3 - uygradn(k,j,1)=(uy(k,i)-uyt(k,i))/delta - uygradn(k,j,2)=(uy(k,i-1)-uyt(k,i-1))/delta - uzgradn(k,j,1)=(uz(k,i)-uzt(k,i))/delta - uzgradn(k,j,2)=(uz(k,i-1)-uzt(k,i-1))/delta - enddo -c write (iout,'(i5,3f8.5,3x,3f8.5,5x,3f8.5,3x,3f8.5)') -c & j,(uzgradt(k,j,1,i),k=1,3),(uzgradn(k,j,1),k=1,3), -c & (uzgradt(k,j,2,i-1),k=1,3),(uzgradn(k,j,2),k=1,3) - enddo - do k=1,3 - dc_norm(k,i)=erij(k) - enddo -cd do k=1,3 -cd write (iout,'(i5,3f8.5,3x,3f8.5,5x,3f8.5,3x,3f8.5)') -cd & k,(uygradt(k,l,1,i),l=1,3),(uygradn(k,l,1),l=1,3), -cd & (uygradt(k,l,2,i-1),l=1,3),(uygradn(k,l,2),l=1,3) -cd write (iout,'(i5,3f8.5,3x,3f8.5,5x,3f8.5,3x,3f8.5)') -cd & k,(uzgradt(k,l,1,i),l=1,3),(uzgradn(k,l,1),l=1,3), -cd & (uzgradt(k,l,2,i-1),l=1,3),(uzgradn(k,l,2),l=1,3) -cd write (iout,'(a)') -cd enddo - enddo - return - end C-------------------------------------------------------------------------- subroutine set_matrices implicit real*8 (a-h,o-z) @@ -2865,7 +3000,7 @@ C-------------------------------------------------------------------------- include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' - include 'COMMON.CONTACTS' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' @@ -2881,18 +3016,26 @@ c write(iout,*) "itype2loc",itype2loc #else do i=3,nres+1 #endif - if (i.gt. nnt+2 .and. i.lt.nct+2) then + ii=ireschain(i-2) +c write (iout,*) "i",i,i-2," ii",ii + if (ii.eq.0) cycle + innt=chain_border(1,ii) + inct=chain_border(2,ii) +c write (iout,*) "i",i,i-2," ii",ii," innt",innt," inct",inct +c if (i.gt. nnt+2 .and. i.lt.nct+2) then + if (i.gt. innt+2 .and. i.lt.inct+2) then iti = itype2loc(itype(i-2)) else iti=nloctyp endif c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then - if (i.gt. nnt+1 .and. i.lt.nct+1) then + if (i.gt. innt+1 .and. i.lt.inct+1) then iti1 = itype2loc(itype(i-1)) else iti1=nloctyp endif -c write(iout,*),i +c write(iout,*),"i",i,i-2," iti",itype(i-2),iti, +c & " iti1",itype(i-1),iti1 #ifdef NEWCORR cost1=dcos(theta(i-1)) sint1=dsin(theta(i-1)) @@ -2958,7 +3101,8 @@ c b2tilde(2,i-2)=-b2(2,i-2) write (iout,*) 'theta=', theta(i-1) #endif #else - if (i.gt. nnt+2 .and. i.lt.nct+2) then + if (i.gt. innt+2 .and. i.lt.inct+2) then +c if (i.gt. nnt+2 .and. i.lt.nct+2) then iti = itype2loc(itype(i-2)) else iti=nloctyp @@ -3003,12 +3147,14 @@ c write(iout,*) 'b2=',(b2(k,i-2),k=1,2) #endif enddo + mu=0.0d0 #ifdef PARMAT do i=ivec_start+2,ivec_end+2 #else do i=3,nres+1 #endif - if (i .lt. nres+1) then +c if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle + if (i .lt. nres+1 .and. itype(i-1).lt.ntyp1) then sin1=dsin(phi(i)) cos1=dcos(phi(i)) sintab(i-2)=sin1 @@ -3045,7 +3191,7 @@ c Ug2(2,1,i-2)=0.0d0 Ug2(2,2,i-2)=0.0d0 endif - if (i .gt. 3 .and. i .lt. nres+1) then + if (i .gt. 3) then obrot_der(1,i-2)=-sin1 obrot_der(2,i-2)= cos1 Ugder(1,1,i-2)= sin1 @@ -3075,7 +3221,8 @@ c Ug2der(2,2,i-2)=0.0d0 endif c if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then - if (i.gt. nnt+2 .and. i.lt.nct+2) then +c if (i.gt. nnt+2 .and. i.lt.nct+2) then + if (i.gt.nnt+2 .and.i.lt.nct+2) then iti = itype2loc(itype(i-2)) else iti=nloctyp @@ -3104,6 +3251,7 @@ c & EE(1,2,iti),EE(2,2,i) c write(iout,*) "Macierz EUG", c & eug(1,1,i-2),eug(1,2,i-2),eug(2,1,i-2), c & eug(2,2,i-2) +#ifdef FOURBODY if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) & then call matmat2(CC(1,1,i-2),Ug(1,1,i-2),CUg(1,1,i-2)) @@ -3112,6 +3260,7 @@ c & eug(2,2,i-2) call matvec2(Ctilde(1,1,i-1),obrot(1,i-2),Ctobr(1,i-2)) call matvec2(Dtilde(1,1,i-2),obrot2(1,i-2),Dtobr2(1,i-2)) endif +#endif else do k=1,2 Ub2(k,i-2)=0.0d0 @@ -3156,6 +3305,7 @@ c mu(k,i-2)=Ub2(k,i-2) cd write (iout,*) 'mu1',mu1(:,i-2) cd write (iout,*) 'mu2',mu2(:,i-2) cd write (iout,*) 'mu',i-2,mu(:,i-2) +#ifdef FOURBODY if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) & then call matmat2(CC(1,1,i-1),Ugder(1,1,i-2),CUgder(1,1,i-2)) @@ -3174,7 +3324,9 @@ C Vectors and matrices dependent on a single virtual-bond dihedral. call matmat2(EUg(1,1,i-2),DD(1,1,i-1),EUgD(1,1,i-2)) call matmat2(EUgder(1,1,i-2),DD(1,1,i-1),EUgDder(1,1,i-2)) endif +#endif enddo +#ifdef FOURBODY C Matrices dependent on two consecutive virtual-bond dihedrals. C The order of matrices is from left to right. if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) @@ -3191,6 +3343,7 @@ c do i=max0(ivec_start,2),ivec_end call matmat2(auxmat(1,1),EUg(1,1,i),Ug2DtEUgder(1,1,1,i)) enddo endif +#endif #if defined(MPI) && defined(PARMAT) #ifdef DEBUG c if (fg_rank.eq.0) then @@ -3259,6 +3412,7 @@ c & " ivec_count",(ivec_count(i),i=0,nfgtasks-1) call MPI_Allgatherv(sintab2(ivec_start),ivec_count(fg_rank1), & MPI_DOUBLE_PRECISION,sintab2(1),ivec_count(0),ivec_displ(0), & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) +#ifdef FOURBODY if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) & then call MPI_Allgatherv(Ctobr(1,ivec_start),ivec_count(fg_rank1), @@ -3334,6 +3488,7 @@ c & " ivec_count",(ivec_count(i),i=0,nfgtasks-1) & MPI_MAT2,Ug2DtEUgder(1,1,1,1),ivec_count(0),ivec_displ(0), & MPI_MAT2,FG_COMM1,IERR) endif +#endif #else c Passes matrix info through the ring isend=fg_rank1 @@ -3378,6 +3533,7 @@ c call flush(iout) & iprev,6600+irecv,FG_COMM,status,IERR) c write (iout,*) "Gather PRECOMP12" c call flush(iout) +#ifdef FOURBODY if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) & then call MPI_SENDRECV(ug2db1t(1,ivec_displ(isend)+1),1, @@ -3397,6 +3553,7 @@ c call flush(iout) & Ug2DtEUgder(1,1,1,ivec_displ(irecv)+1),1, & MPI_PRECOMP23(lenrecv), & iprev,9900+irecv,FG_COMM,status,IERR) +#endif c write (iout,*) "Gather PRECOMP23" c call flush(iout) endif @@ -3449,7 +3606,7 @@ cd enddo cd enddo return end -C-------------------------------------------------------------------------- +C----------------------------------------------------------------------------- subroutine eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) C C This subroutine calculates the average interaction energy and its gradient @@ -3472,7 +3629,11 @@ C include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' +#ifdef FOURBODY include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' +#endif + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' @@ -3545,9 +3706,11 @@ cd enddo eello_turn3=0.0d0 eello_turn4=0.0d0 ind=0 +#ifdef FOURBODY do i=1,nres num_cont_hb(i)=0 enddo +#endif cd print '(a)','Enter EELEC' cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e do i=1,nres @@ -3597,7 +3760,9 @@ c end if num_conti=0 call eelecij(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) +#ifdef FOURBODY num_cont_hb(i)=num_conti +#endif enddo do i=iturn4_start,iturn4_end if (i.lt.1) cycle @@ -3653,12 +3818,16 @@ c endif zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize +#ifdef FOURBODY num_conti=num_cont_hb(i) +#endif c write(iout,*) "JESTEM W PETLI" call eelecij(i,i+3,ees,evdw1,eel_loc) if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) & call eturn4(i,eello_turn4) +#ifdef FOURBODY num_cont_hb(i)=num_conti +#endif enddo ! i C Loop over all neighbouring boxes C do xshift=-1,1 @@ -3725,7 +3894,9 @@ c go to 166 c endif c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) +#ifdef FOURBODY num_conti=num_cont_hb(i) +#endif C I TU KURWA do j=ielstart(i),ielend(i) C do j=16,17 @@ -3741,7 +3912,9 @@ c & .or.itype(j-1).eq.ntyp1 &) cycle call eelecij(i,j,ees,evdw1,eel_loc) enddo ! j +#ifdef FOURBODY num_cont_hb(i)=num_conti +#endif enddo ! i C enddo ! zshift C enddo ! yshift @@ -3759,7 +3932,7 @@ cd print *,"Processor",fg_rank," t_eelecij",t_eelecij end C------------------------------------------------------------------------------- subroutine eelecij(i,j,ees,evdw1,eel_loc) - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifdef MPI include "mpif.h" @@ -3772,21 +3945,44 @@ C------------------------------------------------------------------------------- include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' +#ifdef FOURBODY include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' +#endif + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.TIME1' include 'COMMON.SPLITELE' include 'COMMON.SHIELD' - dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), + double precision 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), & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4), & gmuij2(4),gmuji2(4) + double precision dxi,dyi,dzi + double precision dx_normi,dy_normi,dz_normi,aux + integer j1,j2,lll,num_conti common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + integer k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap,ilist,iresshield + double precision ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp + double precision ees,evdw1,eel_loc,aaa,bbb,ael3i + double precision dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj, + & rij,r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4, + & evdwij,el1,el2,eesij,ees0ij,facvdw,facel,fac1,ecosa, + & ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der, + & a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij, + & ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp, + & ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp, + & ecosgp,ecosam,ecosbm,ecosgm,ghalf,rlocshield + double precision a22,a23,a32,a33,geel_loc_ij,geel_loc_ji + double precision dist_init,xj_safe,yj_safe,zj_safe, + & xj_temp,yj_temp,zj_temp,dist_temp,xmedi,ymedi,zmedi + double precision sscale,sscagrad,scalar + c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ @@ -3890,8 +4086,9 @@ C yj=yj-ymedi C zj=zj-zmedi rij=xj*xj+yj*yj+zj*zj - sss=sscale(sqrt(rij)) - sssgrad=sscagrad(sqrt(rij)) + sss=sscale(sqrt(rij),r_cut_int) + if (sss.eq.0.0d0) return + sssgrad=sscagrad(sqrt(rij),r_cut_int) c if (sss.gt.0.0d0) then rrmij=1.0D0/rij rij=dsqrt(rij) @@ -3926,7 +4123,7 @@ C fac_shield(j)=0.6 fac_shield(i)=1.0 fac_shield(j)=1.0 eesij=(el1+el2) - ees=ees+eesij + ees=ees+eesij*sss endif evdw1=evdw1+evdwij*sss cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') @@ -3935,11 +4132,10 @@ cd & 1.0D0/dsqrt(rrmij),evdwij,eesij, cd & xmedi,ymedi,zmedi,xj,yj,zj if (energy_dec) then - write (iout,'(a6,2i5,0pf7.3,2i5,3e11.3)') - &'evdw1',i,j,evdwij - &,iteli,itelj,aaa,evdw1,sss - write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij, - &fac_shield(i),fac_shield(j) + write (iout,'(a6,2i5,0pf7.3,2i5,e11.3,3f10.5)') + & 'evdw1',i,j,evdwij,iteli,itelj,aaa,evdw1,sss,rij + write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij, + & fac_shield(i),fac_shield(j) endif C @@ -3956,6 +4152,7 @@ C * * Radial derivatives. First process both termini of the fragment (i,j) * + aux=facel+sssgrad*eesij ggg(1)=facel*xj ggg(2)=facel*yj ggg(3)=facel*zj @@ -3992,10 +4189,10 @@ C endif iresshield=shield_list(ilist,j) do k=1,3 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j) - & *2.0 + & *2.0*sss gshieldx(k,iresshield)=gshieldx(k,iresshield)+ & rlocshield - & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0 + & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0*sss gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) @@ -4018,13 +4215,13 @@ C endif do k=1,3 gshieldc(k,i)=gshieldc(k,i)+ - & grad_shield(k,i)*eesij/fac_shield(i)*2.0 + & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss gshieldc(k,j)=gshieldc(k,j)+ - & grad_shield(k,j)*eesij/fac_shield(j)*2.0 + & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss gshieldc(k,i-1)=gshieldc(k,i-1)+ - & grad_shield(k,i)*eesij/fac_shield(i)*2.0 + & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss gshieldc(k,j-1)=gshieldc(k,j-1)+ - & grad_shield(k,j)*eesij/fac_shield(j)*2.0 + & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss enddo endif @@ -4055,15 +4252,10 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo - if (sss.gt.0.0) then - ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj - ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj - ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj - else - ggg(1)=0.0 - ggg(2)=0.0 - ggg(3)=0.0 - endif + facvdw=facvdw+sssgrad*rmij*evdwij + ggg(1)=facvdw*xj + ggg(2)=facvdw*yj + ggg(3)=facvdw*zj c do k=1,3 c ghalf=0.5D0*ggg(k) c gvdwpp(k,i)=gvdwpp(k,i)+ghalf @@ -4084,10 +4276,11 @@ cgrad enddo cgrad enddo #else C MARYSIA - facvdw=(ev1+evdwij)*sss + facvdw=(ev1+evdwij) facel=(el1+eesij) fac1=fac - fac=-3*rrmij*(facvdw+facvdw+facel) + fac=-3*rrmij*(facvdw+facvdw+facel)*sss + & +(evdwij+eesij)*sssgrad*rrmij erij(1)=xj*rmij erij(2)=yj*rmij erij(3)=zj*rmij @@ -4163,11 +4356,11 @@ C print *,"before22", gelc_long(1,i), gelc_long(1,j) 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)) + & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss & *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)) + & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss & *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) @@ -4403,7 +4596,7 @@ 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) + & *fac_shield(i)*fac_shield(j)*sss c if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') c & 'eelloc',i,j,eel_loc_ij C Now derivative over eel_loc @@ -4415,11 +4608,11 @@ C print *,i,j iresshield=shield_list(ilist,i) do k=1,3 rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij - & /fac_shield(i) + & /fac_shield(i)*sss C & *2.0 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+ & rlocshield - & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i) + & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)*sss gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1) & +rlocshield enddo @@ -4428,11 +4621,11 @@ C & *2.0 iresshield=shield_list(ilist,j) do k=1,3 rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij - & /fac_shield(j) + & /fac_shield(j)*sss C & *2.0 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+ & rlocshield - & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j) + & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)*sss gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1) & +rlocshield @@ -4441,13 +4634,13 @@ C & *2.0 do k=1,3 gshieldc_ll(k,i)=gshieldc_ll(k,i)+ - & grad_shield(k,i)*eel_loc_ij/fac_shield(i) + & grad_shield(k,i)*eel_loc_ij/fac_shield(i)*sss gshieldc_ll(k,j)=gshieldc_ll(k,j)+ - & grad_shield(k,j)*eel_loc_ij/fac_shield(j) + & grad_shield(k,j)*eel_loc_ij/fac_shield(j)*sss gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+ - & grad_shield(k,i)*eel_loc_ij/fac_shield(i) + & grad_shield(k,i)*eel_loc_ij/fac_shield(i)*sss gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+ - & grad_shield(k,j)*eel_loc_ij/fac_shield(j) + & grad_shield(k,j)*eel_loc_ij/fac_shield(j)*sss enddo endif @@ -4461,7 +4654,7 @@ C Calculate patrial derivative for theta angle & +a23*gmuij1(2) & +a32*gmuij1(3) & +a33*gmuij1(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss c write(iout,*) "derivative over thatai" c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), c & a33*gmuij1(4) @@ -4477,7 +4670,7 @@ c & a33*gmuij2(4) & +a33*gmuij2(4) gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & geel_loc_ij*wel_loc - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss c Derivative over j residue geel_loc_ji=a22*gmuji1(1) @@ -4490,7 +4683,7 @@ c & a33*gmuji1(4) gloc(nphi+j,icg)=gloc(nphi+j,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss geel_loc_ji= & +a22*gmuji2(1) @@ -4502,7 +4695,7 @@ c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), c & a33*gmuji2(4) gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss #endif cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij @@ -4518,17 +4711,17 @@ C Partial derivatives in virtual-bond dihedral angles gamma & 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) + & *fac_shield(i)*fac_shield(j)*sss 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) + & *fac_shield(i)*fac_shield(j)*sss 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) + & *fac_shield(i)*fac_shield(j)*sss 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) @@ -4544,24 +4737,25 @@ 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) + & *fac_shield(i)*fac_shield(j)*sss 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) + & *fac_shield(i)*fac_shield(j)*sss 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) + & *fac_shield(i)*fac_shield(j)*sss 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) + & *fac_shield(i)*fac_shield(j)*sss enddo ENDIF C Change 12/26/95 to calculate four-body contributions to H-bonding energy c if (j.gt.i+1 .and. num_conti.le.maxconts) then +#ifdef FOURBODY if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0 & .and. num_conti.le.maxconts) then c write (iout,*) i,j," entered corr" @@ -4645,9 +4839,9 @@ 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) + & *fac_shield(i)*fac_shield(j)*sss ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss 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 @@ -4696,11 +4890,17 @@ cd fprimcont=0.0D0 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k) enddo gggp(1)=gggp(1)+ees0pijp*xj + & +ees0p(num_conti,i)/sss*rmij*xj*sssgrad gggp(2)=gggp(2)+ees0pijp*yj + & +ees0p(num_conti,i)/sss*rmij*yj*sssgrad gggp(3)=gggp(3)+ees0pijp*zj + & +ees0p(num_conti,i)/sss*rmij*zj*sssgrad gggm(1)=gggm(1)+ees0mijp*xj + & +ees0m(num_conti,i)/sss*rmij*xj*sssgrad gggm(2)=gggm(2)+ees0mijp*yj + & +ees0m(num_conti,i)/sss*rmij*yj*sssgrad gggm(3)=gggm(3)+ees0mijp*zj + & +ees0m(num_conti,i)/sss*rmij*zj*sssgrad C Derivatives due to the contact function gacont_hbr(1,num_conti,i)=fprimcont*xj gacont_hbr(2,num_conti,i)=fprimcont*yj @@ -4715,28 +4915,28 @@ 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) + & *fac_shield(i)*fac_shield(j)*sss 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) + & *fac_shield(i)*fac_shield(j)*sss gacontp_hb3(k,num_conti,i)=gggp(k) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss 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) + & *fac_shield(i)*fac_shield(j)*sss 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) + & *fac_shield(i)*fac_shield(j)*sss gacontm_hb3(k,num_conti,i)=gggm(k) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss enddo C Diagnostics. Comment out or remove after debugging! @@ -4752,6 +4952,7 @@ cdiag enddo endif ! num_conti.le.maxconts endif ! fcont.gt.0 endif ! j.gt.i+1 +#endif if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then do k=1,4 do l=1,3 @@ -4784,7 +4985,7 @@ C Third- and fourth-order contributions from turns include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' - include 'COMMON.CONTACTS' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' @@ -4967,7 +5168,7 @@ C Third- and fourth-order contributions from turns include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' - include 'COMMON.CONTACTS' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' @@ -5516,7 +5717,7 @@ C This subroutine calculates the excluded-volume interaction energy between C peptide-group centers and side chains and its gradient in virtual-bond and C side-chain vectors. C - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -5529,7 +5730,14 @@ C include 'COMMON.CONTROL' include 'COMMON.SPLITELE' integer xshift,yshift,zshift - dimension ggg(3) + double precision ggg(3) + integer i,iint,j,k,iteli,itypj,subchap + double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1, + & fac,e1,e2,rij + double precision evdw2,evdw2_14,evdwij + double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp, + & dist_temp, dist_init + double precision sscale,sscagrad evdw2=0.0D0 evdw2_14=0.0d0 c print *,boxxsize,boxysize,boxzsize,'wymiary pudla' @@ -5538,7 +5746,7 @@ cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e C do xshift=-1,1 C do yshift=-1,1 C do zshift=-1,1 - if (energy_dec) write (iout,*) "escp:",r_cut,rlamb + if (energy_dec) write (iout,*) "escp:",r_cut_int,rlamb do i=iatscp_s,iatscp_e if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) @@ -5660,11 +5868,11 @@ CHERE IS THE CALCULATION WHICH MIRROR IMAGE IS THE CLOSEST ONE c print *,xj,yj,zj,'polozenie j' rrij=1.0D0/(xj*xj+yj*yj+zj*zj) c print *,rrij - sss=sscale(1.0d0/(dsqrt(rrij))) + sss=sscale(1.0d0/(dsqrt(rrij)),r_cut_int) c print *,r_cut,1.0d0/dsqrt(rrij),sss,'tu patrz' c if (sss.eq.0) print *,'czasem jest OK' if (sss.le.0.0d0) cycle - sssgrad=sscagrad(1.0d0/(dsqrt(rrij))) + sssgrad=sscagrad(1.0d0/(dsqrt(rrij)),r_cut_int) fac=rrij**expon2 e1=fac*fac*aad(itypj,iteli) e2=fac*bad(itypj,iteli) @@ -5675,8 +5883,9 @@ c if (sss.eq.0) print *,'czasem jest OK' endif evdwij=e1+e2 evdw2=evdw2+evdwij*sss - if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)') - & 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli), + if (energy_dec) write (iout,'(a6,2i5,3f7.3,2i3,3e11.3)') + & 'evdw2',i,j,1.0d0/dsqrt(rrij),sss, + & evdwij,iteli,itypj,fac,aad(itypj,iteli), & bad(itypj,iteli) C C Calculate contributions to the gradient in the virtual-bond and SC vectors. @@ -6091,6 +6300,12 @@ c estr=0.0d0 estr1=0.0d0 do i=ibondp_start,ibondp_end +c 3/4/2020 Adam: removed dummy bond graient if Calpha and SC coords are +c used +#ifdef FIVEDIAG + if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) cycle + diff = vbld(i)-vbldp0 +#else if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle c estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax) c do j=1,3 @@ -6101,15 +6316,16 @@ c if (energy_dec) write(iout,*) c & "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax) c else C Checking if it involves dummy (NH3+ or COO-) group - if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then + if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then C YES vbldpDUM is the equlibrium length of spring for Dummy atom - diff = vbld(i)-vbldpDUM - if (energy_dec) write(iout,*) "dum_bond",i,diff - else -C NO vbldp0 is the equlibrium lenght of spring for peptide group - diff = vbld(i)-vbldp0 - endif - if (energy_dec) write (iout,'(a7,i5,4f7.3)') + diff = vbld(i)-vbldpDUM + if (energy_dec) write(iout,*) "dum_bond",i,diff + else +C NO vbldp0 is the equlibrium length of spring for peptide group + diff = vbld(i)-vbldp0 + endif +#endif + if (energy_dec) write (iout,'(a7,i5,4f7.3)') & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff estr=estr+diff*diff do j=1,3 @@ -7061,7 +7277,8 @@ c & sumene4, c & dscp1,dscp2,sumene c sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1)) escloc = escloc + sumene -c write (2,*) "i",i," escloc",sumene,escloc,it,itype(i) + if (energy_dec) write (2,*) "i",i," itype",itype(i)," it",it, + & " escloc",sumene,escloc,it,itype(i) c & ,zz,xx,yy c#define DEBUG #ifdef DEBUG @@ -7443,6 +7660,20 @@ c------------------------------------------------------------------------------ return end c---------------------------------------------------------------------------- +c LICZENIE WIEZOW Z ROWNANIA ENERGII MODELLERA + subroutine e_modeller(ehomology_constr) + ehomology_constr=0.0d0 + write (iout,*) "!!!!!UWAGA, JESTEM W DZIWNEJ PETLI, TEST!!!!!" + return + end +C !!!!!!!! NIE CZYTANE !!!!!!!!!!! + +c------------------------------------------------------------------------------ + subroutine etor_d(etors_d) + etors_d=0.0d0 + return + end +c---------------------------------------------------------------------------- #else subroutine etor(etors) implicit real*8 (a-h,o-z) @@ -7813,6 +8044,637 @@ c do i=1,ndih_constr return end c---------------------------------------------------------------------------- +c MODELLER restraint function + subroutine e_modeller(ehomology_constr) + implicit none + include 'DIMENSIONS' + + double precision ehomology_constr + integer nnn,i,ii,j,k,ijk,jik,ki,kk,nexl,irec,l + integer katy, odleglosci, test7 + real*8 odleg, odleg2, odleg3, kat, kat2, kat3, gdih(max_template) + real*8 Eval,Erot + real*8 distance(max_template),distancek(max_template), + & min_odl,godl(max_template),dih_diff(max_template) + +c +c FP - 30/10/2014 Temporary specifications for homology restraints +c + double precision utheta_i,gutheta_i,sum_gtheta,sum_sgtheta, + & sgtheta + double precision, dimension (maxres) :: guscdiff,usc_diff + double precision, dimension (max_template) :: + & gtheta,dscdiff,uscdiffk,guscdiff2,guscdiff3, + & theta_diff + double precision sum_godl,sgodl,grad_odl3,ggodl,sum_gdih, + & sum_guscdiff,sum_sgdih,sgdih,grad_dih3,usc_diff_i,dxx,dyy,dzz, + & betai,sum_sgodl,dij + double precision dist,pinorm +c + include 'COMMON.SBRIDGE' + include 'COMMON.CHAIN' + include 'COMMON.GEO' + include 'COMMON.DERIV' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.VAR' + include 'COMMON.IOUNITS' +c include 'COMMON.MD' + include 'COMMON.CONTROL' + include 'COMMON.HOMOLOGY' + include 'COMMON.QRESTR' +c +c From subroutine Econstr_back +c + include 'COMMON.NAMES' + include 'COMMON.TIME1' +c + + + do i=1,max_template + distancek(i)=9999999.9 + enddo + + + odleg=0.0d0 + +c Pseudo-energy and gradient from homology restraints (MODELLER-like +c function) +C AL 5/2/14 - Introduce list of restraints +c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d +#ifdef DEBUG + write(iout,*) "------- dist restrs start -------" +#endif + do ii = link_start_homo,link_end_homo + i = ires_homo(ii) + j = jres_homo(ii) + dij=dist(i,j) +c write (iout,*) "dij(",i,j,") =",dij + nexl=0 + do k=1,constr_homology +c write(iout,*) ii,k,i,j,l_homo(k,ii),dij,odl(k,ii) + if(.not.l_homo(k,ii)) then + nexl=nexl+1 + cycle + endif + distance(k)=odl(k,ii)-dij +c write (iout,*) "distance(",k,") =",distance(k) +c +c For Gaussian-type Urestr +c + distancek(k)=0.5d0*distance(k)**2*sigma_odl(k,ii) ! waga_dist rmvd from Gaussian argument +c write (iout,*) "sigma_odl(",k,ii,") =",sigma_odl(k,ii) +c write (iout,*) "distancek(",k,") =",distancek(k) +c distancek(k)=0.5d0*waga_dist*distance(k)**2*sigma_odl(k,ii) +c +c For Lorentzian-type Urestr +c + if (waga_dist.lt.0.0d0) then + sigma_odlir(k,ii)=dsqrt(1/sigma_odl(k,ii)) + distancek(k)=distance(k)**2/(sigma_odlir(k,ii)* + & (distance(k)**2+sigma_odlir(k,ii)**2)) + endif + enddo + +c min_odl=minval(distancek) + do kk=1,constr_homology + if(l_homo(kk,ii)) then + min_odl=distancek(kk) + exit + endif + enddo + do kk=1,constr_homology + if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) + & min_odl=distancek(kk) + enddo + +c write (iout,* )"min_odl",min_odl +#ifdef DEBUG + write (iout,*) "ij dij",i,j,dij + write (iout,*) "distance",(distance(k),k=1,constr_homology) + write (iout,*) "distancek",(distancek(k),k=1,constr_homology) + write (iout,* )"min_odl",min_odl +#endif +#ifdef OLDRESTR + odleg2=0.0d0 +#else + if (waga_dist.ge.0.0d0) then + odleg2=nexl + else + odleg2=0.0d0 + endif +#endif + do k=1,constr_homology +c Nie wiem po co to liczycie jeszcze raz! +c odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/ +c & (2*(sigma_odl(i,j,k))**2)) + if(.not.l_homo(k,ii)) cycle + if (waga_dist.ge.0.0d0) then +c +c For Gaussian-type Urestr +c + godl(k)=dexp(-distancek(k)+min_odl) + odleg2=odleg2+godl(k) +c +c For Lorentzian-type Urestr +c + else + odleg2=odleg2+distancek(k) + endif + +ccc write(iout,779) i,j,k, "odleg2=",odleg2, "odleg3=", odleg3, +ccc & "dEXP(odleg3)=", dEXP(odleg3),"distance(i,j,k)^2=", +ccc & distance(i,j,k)**2, "dist(i+1,j+1)=", dist(i+1,j+1), +ccc & "sigma_odl(i,j,k)=", sigma_odl(i,j,k) + + enddo +c write (iout,*) "godl",(godl(k),k=1,constr_homology) ! exponents +c write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps +#ifdef DEBUG + write (iout,*) "godl",(godl(k),k=1,constr_homology) ! exponents + write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps +#endif + if (waga_dist.ge.0.0d0) then +c +c For Gaussian-type Urestr +c + odleg=odleg-dLOG(odleg2/constr_homology)+min_odl +c +c For Lorentzian-type Urestr +c + else + odleg=odleg+odleg2/constr_homology + endif +c +c write (iout,*) "odleg",odleg ! sum of -ln-s +c Gradient +c +c For Gaussian-type Urestr +c + if (waga_dist.ge.0.0d0) sum_godl=odleg2 + sum_sgodl=0.0d0 + do k=1,constr_homology +c godl=dexp(((-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2)) +c & *waga_dist)+min_odl +c sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist +c + if(.not.l_homo(k,ii)) cycle + if (waga_dist.ge.0.0d0) then +c For Gaussian-type Urestr +c + sgodl=-godl(k)*distance(k)*sigma_odl(k,ii) ! waga_dist rmvd +c +c For Lorentzian-type Urestr +c + else + sgodl=-2*sigma_odlir(k,ii)*(distance(k)/(distance(k)**2+ + & sigma_odlir(k,ii)**2)**2) + endif + sum_sgodl=sum_sgodl+sgodl + +c sgodl2=sgodl2+sgodl +c write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE1" +c write(iout,*) "constr_homology=",constr_homology +c write(iout,*) i, j, k, "TEST K" + enddo + if (waga_dist.ge.0.0d0) then +c +c For Gaussian-type Urestr +c + grad_odl3=waga_homology(iset)*waga_dist + & *sum_sgodl/(sum_godl*dij) +c +c For Lorentzian-type Urestr +c + else +c Original grad expr modified by analogy w Gaussian-type Urestr grad +c grad_odl3=-waga_homology(iset)*waga_dist*sum_sgodl + grad_odl3=-waga_homology(iset)*waga_dist* + & sum_sgodl/(constr_homology*dij) + endif +c +c grad_odl3=sum_sgodl/(sum_godl*dij) + + +c write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE2" +c write(iout,*) (distance(i,j,k)**2), (2*(sigma_odl(i,j,k))**2), +c & (-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2)) + +ccc write(iout,*) godl, sgodl, grad_odl3 + +c grad_odl=grad_odl+grad_odl3 + + do jik=1,3 + ggodl=grad_odl3*(c(jik,i)-c(jik,j)) +ccc write(iout,*) c(jik,i+1), c(jik,j+1), (c(jik,i+1)-c(jik,j+1)) +ccc write(iout,746) "GRAD_ODL_1", i, j, jik, ggodl, +ccc & ghpbc(jik,i+1), ghpbc(jik,j+1) + ghpbc(jik,i)=ghpbc(jik,i)+ggodl + ghpbc(jik,j)=ghpbc(jik,j)-ggodl +ccc write(iout,746) "GRAD_ODL_2", i, j, jik, ggodl, +ccc & ghpbc(jik,i+1), ghpbc(jik,j+1) +c if (i.eq.25.and.j.eq.27) then +c write(iout,*) "jik",jik,"i",i,"j",j +c write(iout,*) "sum_sgodl",sum_sgodl,"sgodl",sgodl +c write(iout,*) "grad_odl3",grad_odl3 +c write(iout,*) "c(",jik,i,")",c(jik,i),"c(",jik,j,")",c(jik,j) +c write(iout,*) "ggodl",ggodl +c write(iout,*) "ghpbc(",jik,i,")", +c & ghpbc(jik,i),"ghpbc(",jik,j,")", +c & ghpbc(jik,j) +c endif + enddo +ccc write(iout,778)"TEST: odleg2=", odleg2, "DLOG(odleg2)=", +ccc & dLOG(odleg2),"-odleg=", -odleg + + enddo ! ii-loop for dist +#ifdef DEBUG + write(iout,*) "------- dist restrs end -------" +c if (waga_angle.eq.1.0d0 .or. waga_theta.eq.1.0d0 .or. +c & waga_d.eq.1.0d0) call sum_gradient +#endif +c Pseudo-energy and gradient from dihedral-angle restraints from +c homology templates +c write (iout,*) "End of distance loop" +c call flush(iout) + kat=0.0d0 +c write (iout,*) idihconstr_start_homo,idihconstr_end_homo +#ifdef DEBUG + write(iout,*) "------- dih restrs start -------" + do i=idihconstr_start_homo,idihconstr_end_homo + write (iout,*) "gloc_init(",i,icg,")",gloc(i,icg) + enddo +#endif + do i=idihconstr_start_homo,idihconstr_end_homo + kat2=0.0d0 +c betai=beta(i,i+1,i+2,i+3) + betai = phi(i) +c write (iout,*) "betai =",betai + do k=1,constr_homology + dih_diff(k)=pinorm(dih(k,i)-betai) +cd write (iout,'(a8,2i4,2f15.8)') "dih_diff",i,k,dih_diff(k) +cd & ,sigma_dih(k,i) +c if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)= +c & -(6.28318-dih_diff(i,k)) +c if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)= +c & 6.28318+dih_diff(i,k) +#ifdef OLD_DIHED + kat3=-0.5d0*dih_diff(k)**2*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument +#else + kat3=(dcos(dih_diff(k))-1)*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument +#endif +c kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i) + gdih(k)=dexp(kat3) + kat2=kat2+gdih(k) +c write(iout,*) "kat2=", kat2, "exp(kat3)=", exp(kat3) +c write(*,*)"" + enddo +c write (iout,*) "gdih",(gdih(k),k=1,constr_homology) ! exps +c write (iout,*) "i",i," betai",betai," kat2",kat2 ! sum of exps +#ifdef DEBUG + write (iout,*) "i",i," betai",betai," kat2",kat2 + write (iout,*) "gdih",(gdih(k),k=1,constr_homology) +#endif + if (kat2.le.1.0d-14) cycle + kat=kat-dLOG(kat2/constr_homology) +c write (iout,*) "kat",kat ! sum of -ln-s + +ccc write(iout,778)"TEST: kat2=", kat2, "DLOG(kat2)=", +ccc & dLOG(kat2), "-kat=", -kat + +c ---------------------------------------------------------------------- +c Gradient +c ---------------------------------------------------------------------- + + sum_gdih=kat2 + sum_sgdih=0.0d0 + do k=1,constr_homology +#ifdef OLD_DIHED + sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i) ! waga_angle rmvd +#else + sgdih=-gdih(k)*dsin(dih_diff(k))*sigma_dih(k,i) ! waga_angle rmvd +#endif +c sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle + sum_sgdih=sum_sgdih+sgdih + enddo +c grad_dih3=sum_sgdih/sum_gdih + grad_dih3=waga_homology(iset)*waga_angle*sum_sgdih/sum_gdih + +c write(iout,*)i,k,gdih,sgdih,beta(i+1,i+2,i+3,i+4),grad_dih3 +ccc write(iout,747) "GRAD_KAT_1", i, nphi, icg, grad_dih3, +ccc & gloc(nphi+i-3,icg) + gloc(i-3,icg)=gloc(i-3,icg)+grad_dih3 +c if (i.eq.25) then +c write(iout,*) "i",i,"icg",icg,"gloc(",i,icg,")",gloc(i,icg) +c endif +ccc write(iout,747) "GRAD_KAT_2", i, nphi, icg, grad_dih3, +ccc & gloc(nphi+i-3,icg) + + enddo ! i-loop for dih +#ifdef DEBUG + write(iout,*) "------- dih restrs end -------" +#endif + +c Pseudo-energy and gradient for theta angle restraints from +c homology templates +c FP 01/15 - inserted from econstr_local_test.F, loop structure +c adapted + +c +c For constr_homology reference structures (FP) +c +c Uconst_back_tot=0.0d0 + Eval=0.0d0 + Erot=0.0d0 +c Econstr_back legacy + do i=1,nres +c do i=ithet_start,ithet_end + dutheta(i)=0.0d0 +c enddo +c do i=loc_start,loc_end + do j=1,3 + duscdiff(j,i)=0.0d0 + duscdiffx(j,i)=0.0d0 + enddo + enddo +c +c do iref=1,nref +c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end +c write (iout,*) "waga_theta",waga_theta + if (waga_theta.gt.0.0d0) then +#ifdef DEBUG + write (iout,*) "usampl",usampl + write(iout,*) "------- theta restrs start -------" +c do i=ithet_start,ithet_end +c write (iout,*) "gloc_init(",nphi+i,icg,")",gloc(nphi+i,icg) +c enddo +#endif +c write (iout,*) "maxres",maxres,"nres",nres + + do i=ithet_start,ithet_end +c +c do i=1,nfrag_back +c ii = ifrag_back(2,i,iset)-ifrag_back(1,i,iset) +c +c Deviation of theta angles wrt constr_homology ref structures +c + utheta_i=0.0d0 ! argument of Gaussian for single k + gutheta_i=0.0d0 ! Sum of Gaussians over constr_homology ref structures +c do j=ifrag_back(1,i,iset)+2,ifrag_back(2,i,iset) ! original loop +c over residues in a fragment +c write (iout,*) "theta(",i,")=",theta(i) + do k=1,constr_homology +c +c dtheta_i=theta(j)-thetaref(j,iref) +c dtheta_i=thetaref(k,i)-theta(i) ! original form without indexing + theta_diff(k)=thetatpl(k,i)-theta(i) +cd write (iout,'(a8,2i4,2f15.8)') "theta_diff",i,k,theta_diff(k) +cd & ,sigma_theta(k,i) + +c + utheta_i=-0.5d0*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta rmvd from Gaussian argument +c utheta_i=-0.5d0*waga_theta*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta? + gtheta(k)=dexp(utheta_i) ! + min_utheta_i? + gutheta_i=gutheta_i+gtheta(k) ! Sum of Gaussians (pk) +c Gradient for single Gaussian restraint in subr Econstr_back +c dutheta(j-2)=dutheta(j-2)+wfrag_back(1,i,iset)*dtheta_i/(ii-1) +c + enddo +c write (iout,*) "gtheta",(gtheta(k),k=1,constr_homology) ! exps +c write (iout,*) "i",i," gutheta_i",gutheta_i ! sum of exps + +c +c Gradient for multiple Gaussian restraint + sum_gtheta=gutheta_i + sum_sgtheta=0.0d0 + do k=1,constr_homology +c New generalized expr for multiple Gaussian from Econstr_back + sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i) ! waga_theta rmvd +c +c sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i)*waga_theta ! right functional form? + sum_sgtheta=sum_sgtheta+sgtheta ! cum variable + enddo +c Final value of gradient using same var as in Econstr_back + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg) + & +sum_sgtheta/sum_gtheta*waga_theta + & *waga_homology(iset) +c dutheta(i-2)=sum_sgtheta/sum_gtheta*waga_theta +c & *waga_homology(iset) +c dutheta(i)=sum_sgtheta/sum_gtheta +c +c Uconst_back=Uconst_back+waga_theta*utheta(i) ! waga_theta added as weight + Eval=Eval-dLOG(gutheta_i/constr_homology) +c write (iout,*) "utheta(",i,")=",utheta(i) ! -ln of sum of exps +c write (iout,*) "Uconst_back",Uconst_back ! sum of -ln-s +c Uconst_back=Uconst_back+utheta(i) + enddo ! (i-loop for theta) +#ifdef DEBUG + write(iout,*) "------- theta restrs end -------" +#endif + endif +c +c Deviation of local SC geometry +c +c Separation of two i-loops (instructed by AL - 11/3/2014) +c +c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end +c write (iout,*) "waga_d",waga_d + +#ifdef DEBUG + write(iout,*) "------- SC restrs start -------" + write (iout,*) "Initial duscdiff,duscdiffx" + do i=loc_start,loc_end + write (iout,*) i,(duscdiff(jik,i),jik=1,3), + & (duscdiffx(jik,i),jik=1,3) + enddo +#endif + do i=loc_start,loc_end + usc_diff_i=0.0d0 ! argument of Gaussian for single k + guscdiff(i)=0.0d0 ! Sum of Gaussians over constr_homology ref structures +c do j=ifrag_back(1,i,iset)+1,ifrag_back(2,i,iset)-1 ! Econstr_back legacy +c write(iout,*) "xxtab, yytab, zztab" +c write(iout,'(i5,3f8.2)') i,xxtab(i),yytab(i),zztab(i) + do k=1,constr_homology +c + dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str? +c Original sign inverted for calc of gradients (s. Econstr_back) + dyy=-yytpl(k,i)+yytab(i) ! ibid y + dzz=-zztpl(k,i)+zztab(i) ! ibid z +c write(iout,*) "dxx, dyy, dzz" +cd write(iout,'(2i5,4f8.2)') k,i,dxx,dyy,dzz,sigma_d(k,i) +c + usc_diff_i=-0.5d0*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d rmvd from Gaussian argument +c usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d? +c uscdiffk(k)=usc_diff(i) + guscdiff2(k)=dexp(usc_diff_i) ! without min_scdiff +c write(iout,*) "i",i," k",k," sigma_d",sigma_d(k,i), +c & " guscdiff2",guscdiff2(k) + guscdiff(i)=guscdiff(i)+guscdiff2(k) !Sum of Gaussians (pk) +c write (iout,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j), +c & xxref(j),yyref(j),zzref(j) + enddo +c +c Gradient +c +c Generalized expression for multiple Gaussian acc to that for a single +c Gaussian in Econstr_back as instructed by AL (FP - 03/11/2014) +c +c Original implementation +c sum_guscdiff=guscdiff(i) +c +c sum_sguscdiff=0.0d0 +c do k=1,constr_homology +c sguscdiff=-guscdiff2(k)*dscdiff(k)*sigma_d(k,i)*waga_d !waga_d? +c sguscdiff=-guscdiff3(k)*dscdiff(k)*sigma_d(k,i)*waga_d ! w min_uscdiff +c sum_sguscdiff=sum_sguscdiff+sguscdiff +c enddo +c +c Implementation of new expressions for gradient (Jan. 2015) +c +c grad_uscdiff=sum_sguscdiff/(sum_guscdiff*dtab) !? + do k=1,constr_homology +c +c New calculation of dxx, dyy, and dzz corrected by AL (07/11), was missing and wrong +c before. Now the drivatives should be correct +c + dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str? +c Original sign inverted for calc of gradients (s. Econstr_back) + dyy=-yytpl(k,i)+yytab(i) ! ibid y + dzz=-zztpl(k,i)+zztab(i) ! ibid z +c +c New implementation +c + sum_guscdiff=guscdiff2(k)*!(dsqrt(dxx*dxx+dyy*dyy+dzz*dzz))* -> wrong! + & sigma_d(k,i) ! for the grad wrt r' +c sum_sguscdiff=sum_sguscdiff+sum_guscdiff +c +c +c New implementation + sum_guscdiff = waga_homology(iset)*waga_d*sum_guscdiff + do jik=1,3 + duscdiff(jik,i-1)=duscdiff(jik,i-1)+ + & sum_guscdiff*(dXX_C1tab(jik,i)*dxx+ + & dYY_C1tab(jik,i)*dyy+dZZ_C1tab(jik,i)*dzz)/guscdiff(i) + duscdiff(jik,i)=duscdiff(jik,i)+ + & sum_guscdiff*(dXX_Ctab(jik,i)*dxx+ + & dYY_Ctab(jik,i)*dyy+dZZ_Ctab(jik,i)*dzz)/guscdiff(i) + duscdiffx(jik,i)=duscdiffx(jik,i)+ + & sum_guscdiff*(dXX_XYZtab(jik,i)*dxx+ + & dYY_XYZtab(jik,i)*dyy+dZZ_XYZtab(jik,i)*dzz)/guscdiff(i) +c +#ifdef DEBUG + write(iout,*) "jik",jik,"i",i + write(iout,*) "dxx, dyy, dzz" + write(iout,'(2i5,3f8.2)') k,i,dxx,dyy,dzz + write(iout,*) "guscdiff2(",k,")",guscdiff2(k) +c write(iout,*) "sum_sguscdiff",sum_sguscdiff +cc write(iout,*) "dXX_Ctab(",jik,i,")",dXX_Ctab(jik,i) +c write(iout,*) "dYY_Ctab(",jik,i,")",dYY_Ctab(jik,i) +c write(iout,*) "dZZ_Ctab(",jik,i,")",dZZ_Ctab(jik,i) +c write(iout,*) "dXX_C1tab(",jik,i,")",dXX_C1tab(jik,i) +c write(iout,*) "dYY_C1tab(",jik,i,")",dYY_C1tab(jik,i) +c write(iout,*) "dZZ_C1tab(",jik,i,")",dZZ_C1tab(jik,i) +c write(iout,*) "dXX_XYZtab(",jik,i,")",dXX_XYZtab(jik,i) +c write(iout,*) "dYY_XYZtab(",jik,i,")",dYY_XYZtab(jik,i) +c write(iout,*) "dZZ_XYZtab(",jik,i,")",dZZ_XYZtab(jik,i) +c write(iout,*) "duscdiff(",jik,i-1,")",duscdiff(jik,i-1) +c write(iout,*) "duscdiff(",jik,i,")",duscdiff(jik,i) +c write(iout,*) "duscdiffx(",jik,i,")",duscdiffx(jik,i) +c endif +#endif + enddo + enddo +c +c uscdiff(i)=-dLOG(guscdiff(i)/(ii-1)) ! Weighting by (ii-1) required? +c usc_diff(i)=-dLOG(guscdiff(i)/constr_homology) ! + min_uscdiff ? +c +c write (iout,*) i," uscdiff",uscdiff(i) +c +c Put together deviations from local geometry + +c Uconst_back=Uconst_back+wfrag_back(1,i,iset)*utheta(i)+ +c & wfrag_back(3,i,iset)*uscdiff(i) + Erot=Erot-dLOG(guscdiff(i)/constr_homology) +c write (iout,*) "usc_diff(",i,")=",usc_diff(i) ! -ln of sum of exps +c write (iout,*) "Uconst_back",Uconst_back ! cum sum of -ln-s +c Uconst_back=Uconst_back+usc_diff(i) +c +c Gradient of multiple Gaussian restraint (FP - 04/11/2014 - right?) +c +c New implment: multiplied by sum_sguscdiff +c + + enddo ! (i-loop for dscdiff) + +c endif + +#ifdef DEBUG + write(iout,*) "------- SC restrs end -------" + write (iout,*) "------ After SC loop in e_modeller ------" + do i=loc_start,loc_end + write (iout,*) "i",i," gradc",(gradc(j,i,icg),j=1,3) + write (iout,*) "i",i," gradx",(gradx(j,i,icg),j=1,3) + enddo + if (waga_theta.eq.1.0d0) then + write (iout,*) "in e_modeller after SC restr end: dutheta" + do i=ithet_start,ithet_end + write (iout,*) i,dutheta(i) + enddo + endif + if (waga_d.eq.1.0d0) then + write (iout,*) "e_modeller after SC loop: duscdiff/x" + do i=1,nres + write (iout,*) i,(duscdiff(j,i),j=1,3) + write (iout,*) i,(duscdiffx(j,i),j=1,3) + enddo + endif +#endif + +c Total energy from homology restraints +#ifdef DEBUG + write (iout,*) "odleg",odleg," kat",kat +#endif +c +c Addition of energy of theta angle and SC local geom over constr_homologs ref strs +c +c ehomology_constr=odleg+kat +c +c For Lorentzian-type Urestr +c + + if (waga_dist.ge.0.0d0) then +c +c For Gaussian-type Urestr +c + ehomology_constr=(waga_dist*odleg+waga_angle*kat+ + & waga_theta*Eval+waga_d*Erot)*waga_homology(iset) +c write (iout,*) "ehomology_constr=",ehomology_constr + else +c +c For Lorentzian-type Urestr +c + ehomology_constr=(-waga_dist*odleg+waga_angle*kat+ + & waga_theta*Eval+waga_d*Erot)*waga_homology(iset) +c write (iout,*) "ehomology_constr=",ehomology_constr + endif +#ifdef DEBUG + write (iout,*) "odleg",waga_dist,odleg," kat",waga_angle,kat, + & "Eval",waga_theta,eval, + & "Erot",waga_d,Erot + write (iout,*) "ehomology_constr",ehomology_constr +#endif + return +c +c FP 01/15 end +c + 748 format(a8,f12.3,a6,f12.3,a7,f12.3) + 747 format(a12,i4,i4,i4,f8.3,f8.3) + 746 format(a12,i4,i4,i4,f8.3,f8.3,f8.3) + 778 format(a7,1X,f10.3,1X,a4,1X,f10.3,1X,a5,1X,f10.3) + 779 format(i3,1X,i3,1X,i2,1X,a7,1X,f7.3,1X,a7,1X,f7.3,1X,a13,1X, + & f7.3,1X,a17,1X,f9.3,1X,a10,1X,f8.3,1X,a10,1X,f8.3) + end +c---------------------------------------------------------------------------- C The rigorous attempt to derive energy function subroutine ebend_kcc(etheta) @@ -7981,6 +8843,7 @@ c write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp return end +#ifdef FOURBODY c---------------------------------------------------------------------------- subroutine multibody(ecorr) C This subroutine calculates multi-body contributions to energy following @@ -7993,6 +8856,8 @@ C contribution equal to sqrt(eps(i,j)*eps(i+1,j+1)) is added. include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' double precision gx(3),gx1(3) logical lprn @@ -8047,6 +8912,8 @@ c------------------------------------------------------------------------------ include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.SHIELD' double precision gx(3),gx1(3) logical lprn @@ -8101,6 +8968,8 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.CONTROL' include 'COMMON.LOCAL' double precision gx(3),gx1(3),time00 @@ -8394,6 +9263,8 @@ c------------------------------------------------------------------------------ parameter (max_cont=maxconts) parameter (max_dim=26) include "COMMON.CONTACTS" + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' double precision zapas(max_dim,maxconts,max_fg_procs), & zapas_recv(max_dim,maxconts,max_fg_procs) common /przechowalnia/ zapas @@ -8465,6 +9336,8 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.CHAIN' include 'COMMON.CONTROL' include 'COMMON.SHIELD' @@ -8835,6 +9708,8 @@ c------------------------------------------------------------------------------ parameter (max_cont=maxconts) parameter (max_dim=70) include "COMMON.CONTACTS" + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' double precision zapas(max_dim,maxconts,max_fg_procs), & zapas_recv(max_dim,maxconts,max_fg_procs) common /przechowalnia/ zapas @@ -8888,6 +9763,8 @@ c------------------------------------------------------------------------------ include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.SHIELD' include 'COMMON.CONTROL' double precision gx(3),gx1(3) @@ -9063,6 +9940,8 @@ C--------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -9128,6 +10007,8 @@ C include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -9514,6 +10395,8 @@ C--------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -9635,6 +10518,8 @@ C--------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -10039,6 +10924,8 @@ c-------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -10179,6 +11066,8 @@ c-------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -10283,6 +11172,8 @@ c---------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -10468,6 +11359,8 @@ c---------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -10583,6 +11476,8 @@ c---------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -10827,6 +11722,8 @@ c---------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -11145,8 +12042,8 @@ cd write (2,*) 'ekont',ekont cd write (2,*) 'eel_turn6',ekont*eel_turn6 return end - C----------------------------------------------------------------------------- +#endif double precision function scalar(u,v) !DIR$ INLINEALWAYS scalar #ifndef OSF @@ -12220,8 +13117,18 @@ c---------------------------------------------------------------------------- include 'COMMON.INTERACT' include 'COMMON.VAR' include 'COMMON.IOUNITS' - include 'COMMON.MD' +c include 'COMMON.MD' +#ifdef LANG0 +#ifdef FIVEDIAG + include 'COMMON.LANGEVIN.lang0.5diag' +#else + include 'COMMON.LANGEVIN.lang0' +#endif +#else + include 'COMMON.LANGEVIN' +#endif include 'COMMON.CONTROL' + include 'COMMON.SAXS' include 'COMMON.NAMES' include 'COMMON.TIME1' include 'COMMON.FFIELD' @@ -12530,8 +13437,18 @@ c---------------------------------------------------------------------------- include 'COMMON.INTERACT' include 'COMMON.VAR' include 'COMMON.IOUNITS' - include 'COMMON.MD' +c include 'COMMON.MD' +#ifdef LANG0 +#ifdef FIVEDIAG + include 'COMMON.LANGEVIN.lang0.5diag' +#else + include 'COMMON.LANGEVIN.lang0' +#endif +#else + include 'COMMON.LANGEVIN' +#endif include 'COMMON.CONTROL' + include 'COMMON.SAXS' include 'COMMON.NAMES' include 'COMMON.TIME1' include 'COMMON.FFIELD'