X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=747b13261294fcdaf24d1f64271bc3ff048ded44;hb=e4035e50115ab9c0e65b0445aed96096a5cb86d5;hp=e6f630d1303b8d4ab6a7acd72c817eda74d05333;hpb=fbad69b401dd3265da402f76b365aa80f6a99135;p=unres.git 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 e6f630d..747b132 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -130,6 +130,8 @@ cmc c if (dyn_ss) call dyn_set_nss c print *,"Processor",myrank," computed USCSC" +c write (iout,*) "SCSC computed OK" +c call flush_(iout) #ifdef TIMING time01=MPI_Wtime() #endif @@ -137,6 +139,14 @@ c print *,"Processor",myrank," computed USCSC" #ifdef TIMING time_vec=time_vec+MPI_Wtime()-time01 #endif +C Introduction of shielding effect first for each peptide group +C the shielding factor is set this factor is describing how each +C peptide group is shielded by side-chains +C the matrix - shield_fac(i) the i index describe the ith between i and i+1 +C write (iout,*) "shield_mode",shield_mode + if (shield_mode.gt.0) then + call set_shield_fac + endif c print *,"Processor",myrank," left VEC_AND_DERIV" if (ipot.lt.6) then #ifdef SPLITELE @@ -163,6 +173,8 @@ c print *,"Processor",myrank," left VEC_AND_DERIV" c call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3, c & eello_turn4) endif +c write (iout,*) "eelec computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed UELEC" C C Calculate excluded-volume interaction energy between peptide groups @@ -179,10 +191,14 @@ C c write (iout,*) "Soft-sphere SCP potential" call escp_soft_sphere(evdw2,evdw2_14) endif +c write (iout,*) "escp computed OK" +c call flush_(iout) c c Calculate the bond-stretching energy c call ebond(estr) +c write (iout,*) "ebond computed OK" +c call flush_(iout) C C Calculate the disulfide-bridge and other energy and the contributions C from other distance constraints. @@ -193,16 +209,21 @@ C C Calculate the virtual-bond-angle energy. C if (wang.gt.0d0) then - call ebend(ebe) + call ebend(ebe,ethetacnstr) else ebe=0 + ethetacnstr=0 endif +c write (iout,*) "ebend computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed UB" C C Calculate the SC local energy. C C print *,"TU DOCHODZE?" call esc(escloc) +c write (iout,*) "esc computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed USC" C C Calculate the virtual-bond torsional energy. @@ -214,6 +235,8 @@ cd print *,'nterm=',nterm etors=0 edihcnstr=0 endif +c write (iout,*) "etor computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed Utor" C C 6/23/01 Calculate double-torsional energy @@ -223,6 +246,8 @@ C else etors_d=0 endif +c write (iout,*) "etor_d computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed Utord" C C 21/5/07 Calculate local sicdechain correlation energy @@ -232,6 +257,8 @@ C else esccor=0.0d0 endif +c write (iout,*) "eback_sc_corr computed OK" +c call flush_(iout) C print *,"PRZED MULIt" c print *,"Processor",myrank," computed Usccorr" C @@ -241,9 +268,13 @@ C n_corr1=0 if ((wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0 & .or. wturn6.gt.0.0d0) .and. ipot.lt.6) then +c write (iout,*) "Calling multibody_eello" +c call flush_(iout) call multibody_eello(ecorr,ecorr5,ecorr6,eturn6,n_corr,n_corr1) -cd write(2,*)'multibody_eello n_corr=',n_corr,' n_corr1=',n_corr1, -cd &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6 +c write(iout,*) +c & 'multibody_eello n_corr=',n_corr,' n_corr1=',n_corr1, +c & " ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6 +c call flush_(iout) else ecorr=0.0d0 ecorr5=0.0d0 @@ -251,9 +282,14 @@ cd &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6 eturn6=0.0d0 endif if ((wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) .and. ipot.lt.6) then +c write (iout,*) "Calling multibody_gb_ecorr" +c call flush_(iout) call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1) -cd write (iout,*) "multibody_hb ecorr",ecorr +c write (iout,*) "Exited multibody_hb ecorr",ecorr +c call flush_(iout) endif +c write (iout,*) "multibody computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed Ucorr" C C If performing constraint dynamics, call the constraint energy @@ -272,7 +308,15 @@ C print *,"przed lipidami" if (wliptran.gt.0) then call Eliptransfer(eliptran) endif -C print *,"za lipidami" +c write (iout,*) "lipid energy computed OK" +c call flush_(iout) + if (AFMlog.gt.0) then + call AFMforce(Eafmforce) + else if (selfguide.gt.0) then + call AFMvel(Eafmforce) + endif +c write (iout,*) "AFMforce computed OK" +c call flush_(iout) #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 #endif @@ -315,11 +359,17 @@ C energia(20)=Uconst+Uconst_back energia(21)=esccor energia(22)=eliptran + energia(23)=Eafmforce + energia(24)=ethetacnstr 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 c print *," Processor",myrank," calls SUM_ENERGY" call sum_energy(energia,.true.) +c write (iout,*) "sum energy OK" +c call flush_(iout) if (dyn_ss) call dyn_set_nss +c write (iout,*) "Exiting energy" +c call flush_(iout) c print *," Processor",myrank," left SUM_ENERGY" #ifdef TIMING time_sumene=time_sumene+MPI_Wtime()-time00 @@ -407,13 +457,16 @@ cMS$ATTRIBUTES C :: proc_proc Uconst=energia(20) esccor=energia(21) eliptran=energia(22) + Eafmforce=energia(23) + ethetacnstr=energia(24) #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d - & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran + & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+Eafmforce + & +ethetacnstr #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc @@ -421,6 +474,8 @@ 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+Uconst+wsccor*esccor+wliptran*eliptran + & +Eafmforce + & +ethetacnstr #endif energia(0)=etot c detecting NaNQ @@ -524,6 +579,8 @@ c enddo & wturn6*gcorr6_turn_long(j,i)+ & wstrain*ghpbc(j,i) & +wliptran*gliptranc(j,i) + & +gradafm(j,i) + & +welec*gshieldc(j,i) enddo enddo @@ -541,6 +598,9 @@ c enddo & wturn6*gcorr6_turn_long(j,i)+ & wstrain*ghpbc(j,i) & +wliptran*gliptranc(j,i) + & +gradafm(j,i) + & +welec*gshieldc(j,i) + enddo enddo #endif @@ -658,6 +718,13 @@ c enddo do i=-1,nct do j=1,3 #ifdef SPLITELE +C print *,gradbufc(1,13) +C print *,welec*gelc(1,13) +C print *,wel_loc*gel_loc(1,13) +C print *,0.5d0*(wscp*gvdwc_scpp(1,13)) +C print *,welec*gelc_long(1,13)+wvdwpp*gvdwpp(1,13) +C print *,wel_loc*gel_loc_long(1,13) +C print *,gradafm(1,13),"AFM" gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ & wel_loc*gel_loc(j,i)+ & 0.5d0*(wscp*gvdwc_scpp(j,i)+ @@ -677,6 +744,11 @@ c enddo & wsccor*gsccorc(j,i) & +wscloc*gscloc(j,i) & +wliptran*gliptranc(j,i) + & +gradafm(j,i) + & +welec*gshieldc(j,i) + & +welec*gshieldc_loc(j,i) + + #else gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ & wel_loc*gel_loc(j,i)+ @@ -697,6 +769,11 @@ c enddo & wsccor*gsccorc(j,i) & +wscloc*gscloc(j,i) & +wliptran*gliptranc(j,i) + & +gradafm(j,i) + & +welec*gshieldc(j,i) + & +welec*gshieldc_loc(j,i) + + #endif gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ @@ -704,6 +781,7 @@ c enddo & wsccor*gsccorx(j,i) & +wscloc*gsclocx(j,i) & +wliptran*gliptranx(j,i) + & +welec*gshieldx(j,i) enddo enddo #ifdef DEBUG @@ -993,15 +1071,17 @@ C------------------------------------------------------------------------ Uconst=energia(20) esccor=energia(21) eliptran=energia(22) + Eafmforce=energia(23) + ethetacnstr=energia(24) #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, & ecorr,wcorr, & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, - & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor, - & edihcnstr,ebr*nss, - & Uconst,eliptran,wliptran,etot + & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr, + & ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, + & etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -1023,10 +1103,13 @@ C------------------------------------------------------------------------ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST= ',1pE16.6,' (Constraint energy)'/ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ + & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ & 'ETOT= ',1pE16.6,' (total)') + #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec, & estr,wbond,ebe,wang, @@ -1034,7 +1117,8 @@ C------------------------------------------------------------------------ & ecorr,wcorr, & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr, - & ebr*nss,Uconst,eliptran,wliptran,etot + & ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, + & etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -1055,9 +1139,11 @@ C------------------------------------------------------------------------ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST=',1pE16.6,' (Constraint energy)'/ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ + & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ & 'ETOT= ',1pE16.6,' (total)') #endif return @@ -1450,6 +1536,7 @@ C include 'COMMON.SBRIDGE' logical lprn integer xshift,yshift,zshift + evdw=0.0D0 ccccc energy_dec=.false. C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -1546,10 +1633,36 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN + +c write(iout,*) "PRZED ZWYKLE", evdwij call dyn_ssbond_ene(i,j,evdwij) +c write(iout,*) "PO ZWYKLE", evdwij + evdw=evdw+evdwij if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') & 'evdw',i,j,evdwij,' ss' +C triple bond artifac removal + do k=j+1,iend(i,iint) +C search over all next residues + if (dyn_ss_mask(k)) then +C check if they are cysteins +C write(iout,*) 'k=',k + +c write(iout,*) "PRZED TRI", evdwij + evdwij_przed_tri=evdwij + call triple_ssbond_ene(i,j,k,evdwij) +c if(evdwij_przed_tri.ne.evdwij) then +c write (iout,*) "TRI:", evdwij, evdwij_przed_tri +c endif + +c write(iout,*) "PO TRI", evdwij +C call the energy function that removes the artifical triple disulfide +C bond the soubroutine is located in ssMD.F + evdw=evdw+evdwij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') + & 'evdw',i,j,evdwij,'tss' + endif!dyn_ss_mask(k) + enddo! k ELSE ind=ind+1 itypj=iabs(itype(j)) @@ -1623,7 +1736,7 @@ C what fraction I am in C lipbufthick is thickenes of lipid buffore sslipj=sscalelip(fracinbuf) ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick - elseif (zi.gt.bufliptop) then + elseif (zj.gt.bufliptop) then fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) sslipj=sscalelip(fracinbuf) ssgradlipj=sscagradlip(fracinbuf)/lipbufthick @@ -1825,12 +1938,12 @@ C the energy transfer exist if (zi.lt.buflipbot) then C what fraction I am in fracinbuf=1.0d0- - & ((positi-bordlipbot)/lipbufthick) + & ((zi-bordlipbot)/lipbufthick) C lipbufthick is thickenes of lipid buffore sslipi=sscalelip(fracinbuf) ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick elseif (zi.gt.bufliptop) then - fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) + fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick) sslipi=sscalelip(fracinbuf) ssgradlipi=sscagradlip(fracinbuf)/lipbufthick else @@ -1893,12 +2006,12 @@ C the energy transfer exist if (zj.lt.buflipbot) then C what fraction I am in fracinbuf=1.0d0- - & ((positi-bordlipbot)/lipbufthick) + & ((zj-bordlipbot)/lipbufthick) C lipbufthick is thickenes of lipid buffore sslipj=sscalelip(fracinbuf) ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick - elseif (zi.gt.bufliptop) then - fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) + elseif (zj.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) sslipj=sscalelip(fracinbuf) ssgradlipj=sscagradlip(fracinbuf)/lipbufthick else @@ -2700,12 +2813,37 @@ c write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2) c write(iout,*) 'b1=',b1(1,i-2) c write (iout,*) 'theta=', theta(i-1) enddo +#else + if (i.gt. nnt+2 .and. i.lt.nct+2) then + iti = itortyp(itype(i-2)) + else + iti=ntortyp+1 + endif +c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then + if (i.gt. nnt+1 .and. i.lt.nct+1) then + iti1 = itortyp(itype(i-1)) + else + iti1=ntortyp+1 + endif + b1(1,i-2)=b(3,iti) + b1(2,i-2)=b(5,iti) + b2(1,i-2)=b(2,iti) + b2(2,i-2)=b(4,iti) + b1tilde(1,i-2)=b1(1,i-2) + b1tilde(2,i-2)=-b1(2,i-2) + b2tilde(1,i-2)=b2(1,i-2) + b2tilde(2,i-2)=-b2(2,i-2) + EE(1,2,i-2)=eeold(1,2,iti) + EE(2,1,i-2)=eeold(2,1,iti) + EE(2,2,i-2)=eeold(2,2,iti) + EE(1,1,i-2)=eeold(1,1,iti) + enddo +#endif #ifdef PARMAT do i=ivec_start+2,ivec_end+2 #else do i=3,nres+1 #endif -#endif if (i .lt. nres+1) then sin1=dsin(phi(i)) cos1=dcos(phi(i)) @@ -2841,6 +2979,7 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then do k=1,2 mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1) enddo +C write (iout,*) 'mumu',i,b1(1,i-1),Ub2(1,i-2) c write (iout,*) 'mu ',mu(:,i-2),i-2 cd write (iout,*) 'mu1',mu1(:,i-2) cd write (iout,*) 'mu2',mu2(:,i-2) @@ -3253,11 +3392,18 @@ C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition if (i.le.1) cycle C write(iout,*) "tu jest i",i if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 +C changes suggested by Ana to avoid out of bounds + & .or.((i+4).gt.nres) + & .or.((i-1).le.0) +C end of changes by Ana & .or. itype(i+2).eq.ntyp1 - & .or. itype(i+3).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 - & .or. itype(i+4).eq.ntyp1 - & ) cycle + & .or. itype(i+3).eq.ntyp1) cycle + if(i.gt.1)then + if(itype(i-1).eq.ntyp1)cycle + end if + if(i.LT.nres-3)then + if (itype(i+4).eq.ntyp1) cycle + end if dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -3281,6 +3427,10 @@ C write(iout,*) "tu jest i",i do i=iturn4_start,iturn4_end if (i.le.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 +C changes suggested by Ana to avoid out of bounds + & .or.((i+5).gt.nres) + & .or.((i-1).le.0) +C end of changes suggested by Ana & .or. itype(i+3).eq.ntyp1 & .or. itype(i+4).eq.ntyp1 & .or. itype(i+5).eq.ntyp1 @@ -3342,9 +3492,15 @@ C do zshift=-1,1 c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c +CTU KURWA do i=iatel_s,iatel_e +C do i=75,75 if (i.le.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 +C changes suggested by Ana to avoid out of bounds + & .or.((i+2).gt.nres) + & .or.((i-1).le.0) +C end of changes by Ana & .or. itype(i+2).eq.ntyp1 & .or. itype(i-1).eq.ntyp1 & ) cycle @@ -3395,10 +3551,16 @@ c endif c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) +C I TU KURWA do j=ielstart(i),ielend(i) +C do j=16,17 C write (iout,*) i,j if (j.le.1) cycle if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1 +C changes suggested by Ana to avoid out of bounds + & .or.((j+2).gt.nres) + & .or.((j-1).le.0) +C end of changes by Ana & .or.itype(j+2).eq.ntyp1 & .or.itype(j-1).eq.ntyp1 &) cycle @@ -3441,6 +3603,7 @@ C------------------------------------------------------------------------------- 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), & 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), @@ -3573,10 +3736,22 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)) el2=fac4*fac C MARYSIA - eesij=(el1+el2) +C eesij=(el1+el2) C 12/26/95 - for the evaluation of multi-body H-bonding interactions ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg) + if (shield_mode.gt.0) then +C fac_shield(i)=0.4 +C fac_shield(j)=0.6 + el1=el1*fac_shield(i)*fac_shield(j) + el2=el2*fac_shield(i)*fac_shield(j) + eesij=(el1+el2) + ees=ees+eesij + else + fac_shield(i)=1.0 + fac_shield(j)=1.0 + eesij=(el1+el2) ees=ees+eesij + endif 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, @@ -3600,22 +3775,99 @@ C erij(1)=xj*rmij erij(2)=yj*rmij erij(3)=zj*rmij + * * Radial derivatives. First process both termini of the fragment (i,j) * 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) + gshieldx(k,iresshield)=gshieldx(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) + gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield +C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield) +C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) +C if (iresshield.gt.i) then +C do ishi=i+1,iresshield-1 +C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield +C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) +C +C enddo +C else +C do ishi=iresshield,i +C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield +C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) +C +C enddo +C endif + enddo + enddo + do ilist=1,ishield_list(j) + iresshield=shield_list(ilist,j) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j) + gshieldx(k,iresshield)=gshieldx(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) + gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield + +C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) +C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield) +C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) +C if (iresshield.gt.j) then +C do ishi=j+1,iresshield-1 +C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield +C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) +C +C enddo +C else +C do ishi=iresshield,j +C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield +C & -grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) +C enddo +C endif + enddo + enddo + + do k=1,3 + gshieldc(k,i)=gshieldc(k,i)+ + & grad_shield(k,i)*eesij/fac_shield(i) + gshieldc(k,j)=gshieldc(k,j)+ + & grad_shield(k,j)*eesij/fac_shield(j) + gshieldc(k,i-1)=gshieldc(k,i-1)+ + & grad_shield(k,i)*eesij/fac_shield(i) + gshieldc(k,j-1)=gshieldc(k,j-1)+ + & grad_shield(k,j)*eesij/fac_shield(j) + + enddo + endif c do k=1,3 c ghalf=0.5D0*ggg(k) c gelc(k,i)=gelc(k,i)+ghalf c gelc(k,j)=gelc(k,j)+ghalf c enddo c 9/28/08 AL Gradient compotents will be summed only at the end +C print *,"before", gelc_long(1,i), gelc_long(1,j) do k=1,3 gelc_long(k,j)=gelc_long(k,j)+ggg(k) +C & +grad_shield(k,j)*eesij/fac_shield(j) gelc_long(k,i)=gelc_long(k,i)-ggg(k) +C & +grad_shield(k,i)*eesij/fac_shield(i) +C gelc_long(k,i-1)=gelc_long(k,i-1) +C & +grad_shield(k,i)*eesij/fac_shield(i) +C gelc_long(k,j-1)=gelc_long(k,j-1) +C & +grad_shield(k,j)*eesij/fac_shield(j) enddo +C print *,"bafter", gelc_long(1,i), gelc_long(1,j) + * * Loop over residues i+1 thru j-1. * @@ -3664,8 +3916,11 @@ C MARYSIA * Radial derivatives. First process both termini of the fragment (i,j) * ggg(1)=fac*xj +C+eesij*grad_shield(1,i)+eesij*grad_shield(1,j) ggg(2)=fac*yj +C+eesij*grad_shield(2,i)+eesij*grad_shield(2,j) ggg(3)=fac*zj +C+eesij*grad_shield(3,i)+eesij*grad_shield(3,j) c do k=1,3 c ghalf=0.5D0*ggg(k) c gelc(k,i)=gelc(k,i)+ghalf @@ -3708,7 +3963,8 @@ c 9/28/08 AL Gradient compotents will be summed only at the end cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3), cd & (dcosg(k),k=1,3) do k=1,3 - ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k) + ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))* + & fac_shield(i)*fac_shield(j) enddo c do k=1,3 c ghalf=0.5D0*ggg(k) @@ -3724,16 +3980,21 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo +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) + & +((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)*fac_shield(j) gelc(k,j)=gelc(k,j) - & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) - & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) + & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) + & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)) + & *fac_shield(i)*fac_shield(j) gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo +C print *,"before33", gelc_long(1,i), gelc_long(1,j) + C MARYSIA c endif !sscale IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 @@ -3936,7 +4197,7 @@ C Contribution to the local-electrostatic energy coming from the i-j pair & +a33*muij(4) c write (iout,*) 'i',i,' j',j,itype(i),itype(j), c & ' eel_loc_ij',eel_loc_ij -c write(iout,*) 'muije=',muij(1),muij(2),muij(3),muij(4) +C write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4) C Calculate patrial derivative for theta angle #ifdef NEWCORR geel_loc_ij=a22*gmuij1(1) @@ -5058,8 +5319,13 @@ C include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' + include 'COMMON.CONTROL' dimension ggg(3) ehpb=0.0D0 + do i=1,3 + ggg(i)=0.0d0 + enddo +C write (iout,*) ,"link_end",link_end,constr_dist cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr cd write(iout,*)'link_start=',link_start,' link_end=',link_end if (link_end.eq.0) return @@ -5086,27 +5352,84 @@ cmc if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds if (.not.dyn_ss .and. i.le.nss) then C 15/02/13 CC dynamic SSbond - additional check - if (ii.gt.nres - & .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then + if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. + & iabs(itype(jjj)).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij endif cd write (iout,*) "eij",eij +cd & ' waga=',waga,' fac=',fac + else if (ii.gt.nres .and. jj.gt.nres) then +c Restraints from contact prediction + dd=dist(ii,jj) + if (constr_dist.eq.11) then + ehpb=ehpb+fordepth(i)**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd + if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, + & ehpb,fordepth(i),dd + else + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "beta nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else + dd=dist(ii,jj) + rdis=dd-dhpb(i) +C Get the force constant corresponding to this distance. + waga=forcon(i) +C Calculate the contribution to energy. + ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "beta reg",dd,waga*rdis*rdis +C +C Evaluate gradient. +C + fac=waga*rdis/dd + endif + endif + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo + do j=1,3 + ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) + ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) + enddo + do k=1,3 + ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) + ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) + enddo else C Calculate the distance between the two points and its difference from the C target distance. dd=dist(ii,jj) + if (constr_dist.eq.11) then + ehpb=ehpb+fordepth(i)**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd + if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, + & ehpb,fordepth(i),dd + else + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "alph nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else rdis=dd-dhpb(i) C Get the force constant corresponding to this distance. waga=forcon(i) C Calculate the contribution to energy. ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "alpha reg",dd,waga*rdis*rdis C C Evaluate gradient. C fac=waga*rdis/dd -cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, -cd & ' waga=',waga,' fac=',fac + endif + endif do j=1,3 ggg(j)=fac*(c(j,jj)-c(j,ii)) enddo @@ -5130,7 +5453,7 @@ cgrad enddo enddo endif enddo - ehpb=0.5D0*ehpb + if (constr_dist.ne.11) ehpb=0.5D0*ehpb return end C-------------------------------------------------------------------------- @@ -5321,7 +5644,7 @@ c end #ifdef CRYST_THETA C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@ -5338,6 +5661,7 @@ C include 'COMMON.NAMES' include 'COMMON.FFIELD' include 'COMMON.CONTROL' + include 'COMMON.TORCNSTR' common /calcthet/ term1,term2,termm,diffak,ratak, & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq, & delthe0,sig0inv,sigtc,sigsqtc,delthec,it @@ -5348,6 +5672,8 @@ c time12=1.0d0 etheta=0.0D0 c write (*,'(a,i2)') 'EBEND ICG=',icg do i=ithet_start,ithet_end +c write (iout,*) "ebend: i=",i +c call flush_(iout) if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 & .or.itype(i).eq.ntyp1) cycle C Zero the energy function and its derivative at 0 or pi. @@ -5449,6 +5775,39 @@ C Derivatives of the "mean" values in gamma1 and gamma2. if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2 gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg) enddo +c write (iout,*) "Exit loop" +c call flush_(iout) + ethetacnstr=0.0d0 +c write (iout,*) ithetaconstr_start,ithetaconstr_end,"TU" +c call flush_(iout) + do i=max0(ithetaconstr_start,1),ithetaconstr_end + itheta=itheta_constr(i) + thetiii=theta(itheta) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetacnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetacnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif + if (energy_dec) then + write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", + & i,itheta,rad2deg*thetiii, + & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), + & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, + & gloc(itheta+nphi-2,icg) + endif + enddo +c write (iout,*) "Exit ebend" +c call flush_(iout) + C Ufff.... We've done all this!!! return end @@ -5565,7 +5924,7 @@ C "Thank you" to MAPLE (probably spared one day of hand-differentiation). end #else C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@ -5584,6 +5943,7 @@ C include 'COMMON.NAMES' include 'COMMON.FFIELD' include 'COMMON.CONTROL' + include 'COMMON.TORCNSTR' double precision coskt(mmaxtheterm),sinkt(mmaxtheterm), & cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle), & sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble), @@ -5591,11 +5951,11 @@ C logical lprn /.false./, lprn1 /.false./ etheta=0.0D0 do i=ithet_start,ithet_end + if (i.le.2) cycle c print *,i,itype(i-1),itype(i),itype(i-2) if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 & .or.itype(i).eq.ntyp1) cycle -C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF - +C print *,i,theta(i) if (iabs(itype(i+1)).eq.20) iblock=2 if (iabs(itype(i+1)).ne.20) iblock=1 dethetai=0.0d0 @@ -5607,6 +5967,16 @@ C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo +C print *,ethetai + if (i.eq.3) then + phii=0.0d0 + ityp1=nthetyp+1 + do k=1,nsingle + cosph1(k)=0.0d0 + sinph1(k)=0.0d0 + enddo + else + if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) @@ -5622,12 +5992,13 @@ C propagation of chirality for glycine type enddo else phii=0.0d0 - ityp1=nthetyp+1 do k=1,nsingle + ityp1=ithetyp((itype(i-2))) cosph1(k)=0.0d0 sinph1(k)=0.0d0 enddo endif + endif if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) @@ -5643,7 +6014,7 @@ C propagation of chirality for glycine type enddo else phii1=0.0d0 - ityp3=nthetyp+1 + ityp3=ithetyp((itype(i))) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 @@ -5693,6 +6064,7 @@ C propagation of chirality for glycine type enddo write(iout,*) "ethetai",ethetai endif +C print *,ethetai do m=1,ntheterm2 do k=1,nsingle aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k) @@ -5713,10 +6085,16 @@ C propagation of chirality for glycine type & ccthet(k,m,ityp1,ityp2,ityp3,iblock)," ddthet", & ddthet(k,m,ityp1,ityp2,ityp3,iblock)," eethet", & eethet(k,m,ityp1,ityp2,ityp3,iblock)," ethetai",ethetai +C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k) enddo enddo +C print *,"cosph1", (cosph1(k), k=1,nsingle) +C print *,"cosph2", (cosph2(k), k=1,nsingle) +C print *,"sinph1", (sinph1(k), k=1,nsingle) +C print *,"sinph2", (sinph2(k), k=1,nsingle) if (lprn) & write(iout,*) "ethetai",ethetai +C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k) do m=1,ntheterm3 do k=2,ndouble do l=1,k-1 @@ -5752,6 +6130,7 @@ C propagation of chirality for glycine type enddo 10 continue c lprn1=.true. +C print *,ethetai if (lprn1) & write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') & i,theta(i)*rad2deg,phii*rad2deg, @@ -5760,8 +6139,37 @@ c lprn1=.false. etheta=etheta+ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 - gloc(nphi+i-2,icg)=wang*dethetai+gloc(nphi+i-2,icg) + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai enddo +C now constrains + ethetacnstr=0.0d0 +C print *,ithetaconstr_start,ithetaconstr_end,"TU" + do i=max0(ithetaconstr_start,1),ithetaconstr_end + itheta=itheta_constr(i) + thetiii=theta(itheta) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif + if (energy_dec) then + write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", + & i,itheta,rad2deg*thetiii, + & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), + & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, + & gloc(itheta+nphi-2,icg) + endif + enddo + return end #endif @@ -6581,12 +6989,12 @@ c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) difi=phii-phi0(i) if (difi.gt.drange(i)) then difi=difi-drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 else if (difi.lt.-drange(i)) then difi=difi+drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)**difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 endif ! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii, ! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) @@ -6619,7 +7027,7 @@ c---------------------------------------------------------------------------- logical lprn C Set lprn=.true. for debugging lprn=.false. -c lprn=.true. +c lprn=.true. etors=0.0D0 do i=iphi_start,iphi_end C ANY TWO ARE DUMMY ATOMS in row CYCLE @@ -6692,18 +7100,21 @@ c do i=1,ndih_constr difi=pinorm(phii-phi0(i)) if (difi.gt.drange(i)) then difi=difi-drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 else if (difi.lt.-drange(i)) then difi=difi+drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 else difi=0.0 endif -cd write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii, -cd & rad2deg*phi0(i), rad2deg*drange(i), -cd & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) + if (energy_dec) then + write (iout,'(a6,2i5,4f8.3,2e14.5)') "edihc", + & i,itori,rad2deg*phii, + & rad2deg*phi0(i), rad2deg*drange(i), + & rad2deg*difi,0.25d0*ftors(i)*difi**4,gloc(itori-3,icg) + endif enddo cd write (iout,*) 'edihcnstr',edihcnstr return @@ -7009,7 +7420,7 @@ C Set lprn=.true. for debugging if (lprn) then write (iout,'(a)') 'Contact function values before RECEIVE:' do i=nnt,nct-2 - write (iout,'(2i3,50(1x,i2,f5.2))') + write (iout,'(2i3,50(1x,i3,f5.2))') & i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), & j=1,num_cont_hb(i)) enddo @@ -7248,6 +7659,7 @@ C Calculate the local-electrostatic correlation terms jp1=iabs(j1) c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, c & ' jj=',jj,' kk=',kk +c call flush(iout) if ((j.gt.0 .and. j1.gt.0 .or. j.gt.0 .and. j1.lt.0 & .or. j.lt.0 .and. j1.gt.0) .and. & (jp1.eq.jp+1 .or. jp1.eq.jp-1)) then @@ -7265,8 +7677,9 @@ c ecorr=ecorr+ehbcorr(i,j,i+1,j,jj,kk,0.60D0,-0.40D0) enddo ! kk do kk=1,num_conti j1=jcont_hb(kk,i) -c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, -c & ' jj=',jj,' kk=',kk +c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, +c & ' jj=',jj,' kk=',kk +c call flush(iout) if (j1.eq.j+1) then C Contacts I-J and (I+1)-J occur simultaneously. C The system loses extra energy. @@ -7360,6 +7773,7 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding include 'COMMON.CONTACTS' include 'COMMON.CHAIN' include 'COMMON.CONTROL' + include 'COMMON.TORSION' double precision gx(3),gx1(3) integer num_cont_hb_old(maxres) logical lprn,ldone @@ -7368,6 +7782,8 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding C Set lprn=.true. for debugging lprn=.false. eturn6=0.0d0 +c write (iout,*) "MULTIBODY_EELLO" +c call flush(iout) #ifdef MPI do i=1,nres num_cont_hb_old(i)=num_cont_hb(i) @@ -7378,12 +7794,12 @@ C Set lprn=.true. for debugging if (lprn) then write (iout,'(a)') 'Contact function values before RECEIVE:' do i=nnt,nct-2 - write (iout,'(2i3,50(1x,i2,f5.2))') + write (iout,'(2i3,50(1x,i3,f5.2))') & i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), & j=1,num_cont_hb(i)) enddo + call flush(iout) endif - call flush(iout) do i=1,ntask_cont_from ncont_recv(i)=0 enddo @@ -7585,10 +8001,15 @@ c call flush(iout) if (lprn) then write (iout,'(a)') 'Contact function values:' do i=nnt,nct-2 - write (iout,'(2i3,50(1x,i2,5f6.3))') + write (iout,'(2i3,50(1x,i3,5f6.3))') & i,num_cont_hb(i),(jcont_hb(j,i),d_cont(j,i), & ((a_chuj(ll,kk,j,i),ll=1,2),kk=1,2),j=1,num_cont_hb(i)) enddo + write (iout,*) "itortyp" + do i=1,nres + write (iout,*) i,itype(i),itortyp(itype(i)) + enddo + call flush(iout) endif ecorr=0.0D0 ecorr5=0.0d0 @@ -7629,8 +8050,11 @@ c write (iout,*) "corr loop i",i do kk=1,num_conti1 j1=jcont_hb(kk,i1) jp1=iabs(j1) -c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, -c & ' jj=',jj,' kk=',kk + if (lprn) then + write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, + & ' jj=',jj,' kk=',kk + call flush(iout) + endif c if (j1.eq.j+1 .or. j1.eq.j-1) then if ((j.gt.0 .and. j1.gt.0 .or. j.gt.0 .and. j1.lt.0 & .or. j.lt.0 .and. j1.gt.0) .and. @@ -7670,6 +8094,12 @@ c do iii=1,nres c write (iout,'(i5,3f10.5)') c & iii,(gradcorr5(jjj,iii),jjj=1,3) c enddo +c write (iout,*) "ecorr4" +c call flush(iout) +c write (iout,*) "eello5:",i,jp,i+1,jp1,jj,kk, +c & itype(jp),itype(i+1),itype(jp1), +c & itortyp(itype(jp)),itortyp(itype(i+1)),itortyp(itype(jp1)) +c call flush(iout) if (wcorr5.gt.0.0d0) & ecorr5=ecorr5+eello5(i,jp,i+1,jp1,jj,kk) c write (iout,*) "gradcorr5 after eello5" @@ -7682,12 +8112,16 @@ c enddo 2 'ecorr5',i,j,i+1,j1,eello5(i,jp,i+1,jp1,jj,kk) cd write(2,*)'wcorr6',wcorr6,' wturn6',wturn6 cd write(2,*)'ijkl',i,jp,i+1,jp1 +c write (iout,*) "ecorr5" +c call flush(iout) if (wcorr6.gt.0.0d0 .and. (jp.ne.i+4 .or. jp1.ne.i+3 & .or. wturn6.eq.0.0d0))then cd write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1 ecorr6=ecorr6+eello6(i,jp,i+1,jp1,jj,kk) if (energy_dec) write (iout,'(a6,4i5,0pf7.3)') 1 'ecorr6',i,j,i+1,j1,eello6(i,jp,i+1,jp1,jj,kk) +c write (iout,*) "ecorr6" +c call flush(iout) cd write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5, cd & 'ecorr6=',ecorr6 cd write (iout,'(4e15.5)') sred_geom, @@ -7701,10 +8135,13 @@ cd write (iout,*) '******eturn6: i,j,i+1,j1',i,jip,i+1,jp1 if (energy_dec) write (iout,'(a6,4i5,0pf7.3)') 1 'eturn6',i,j,i+1,j1,eello_turn6(i,jj,kk) cd write (2,*) 'multibody_eello:eturn6',eturn6 +c write (iout,*) "ecorr4" +c call flush(iout) endif ENDIF 1111 continue endif + if (energy_dec) call flush(iout) enddo ! kk enddo ! jj enddo ! i @@ -8467,9 +8904,9 @@ cd endif cd write (iout,*) cd & 'EELLO5: Contacts have occurred for peptide groups',i,j, cd & ' and',k,l - itk=itortyp(itype(k)) - itl=itortyp(itype(l)) - itj=itortyp(itype(j)) +c itk=itortyp(itype(k)) +c itl=itortyp(itype(l)) +c itj=itortyp(itype(j)) eello5_1=0.0d0 eello5_2=0.0d0 eello5_3=0.0d0 @@ -8538,7 +8975,7 @@ C Cartesian gradient c goto 1112 c1111 continue C Contribution from graph II - call transpose2(EE(1,1,itk),auxmat(1,1)) + call transpose2(EE(1,1,k),auxmat(1,1)) call matmat2(auxmat(1,1),AEA(1,1,1),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -8619,7 +9056,7 @@ C Cartesian gradient cd goto 1112 C Contribution from graph IV cd1110 continue - call transpose2(EE(1,1,itl),auxmat(1,1)) + call transpose2(EE(1,1,l),auxmat(1,1)) call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -8692,7 +9129,7 @@ C Cartesian gradient cd goto 1112 C Contribution from graph IV 1110 continue - call transpose2(EE(1,1,itj),auxmat(1,1)) + call transpose2(EE(1,1,j),auxmat(1,1)) call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -8780,9 +9217,9 @@ cd ghalf=0.0d0 cold ghalf=0.5d0*eel5*eij*gacont_hbr(ll,kk,k) cgrad ghalf=0.5d0*ggg2(ll) cd ghalf=0.0d0 - gradcorr5(ll,k)=gradcorr5(ll,k)+ghalf+ekont*derx(ll,2,2) + gradcorr5(ll,k)=gradcorr5(ll,k)+ekont*derx(ll,2,2) gradcorr5(ll,k+1)=gradcorr5(ll,k+1)+ekont*derx(ll,3,2) - gradcorr5(ll,l)=gradcorr5(ll,l)+ghalf+ekont*derx(ll,4,2) + gradcorr5(ll,l)=gradcorr5(ll,l)+ekont*derx(ll,4,2) gradcorr5(ll,l1)=gradcorr5(ll,l1)+ekont*derx(ll,5,2) gradcorr5_long(ll,l)=gradcorr5_long(ll,l)+gradcorr5kl gradcorr5_long(ll,k)=gradcorr5_long(ll,k)-gradcorr5kl @@ -8989,7 +9426,7 @@ C o o o o C C i i C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC - itk=itortyp(itype(k)) +c itk=itortyp(itype(k)) s1= scalar2(AEAb1(1,2,imat),CUgb2(1,i)) s2=-scalar2(AEAb2(1,1,imat),Ug2Db1t(1,k)) s3= scalar2(AEAb2(1,1,imat),CUgb2(1,k)) @@ -9279,19 +9716,19 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C 4/7/01 AL Component s1 was removed, because it pertains to the respective C energy moment and not to the cluster cumulant. - iti=itortyp(itype(i)) - if (j.lt.nres-1) then - itj1=itortyp(itype(j+1)) - else - itj1=ntortyp - endif - itk=itortyp(itype(k)) - itk1=itortyp(itype(k+1)) - if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) - else - itl1=ntortyp - endif +c iti=itortyp(itype(i)) +c if (j.lt.nres-1) then +c itj1=itortyp(itype(j+1)) +c else +c itj1=ntortyp +c endif +c itk=itortyp(itype(k)) +c itk1=itortyp(itype(k+1)) +c if (l.lt.nres-1) then +c itl1=itortyp(itype(l+1)) +c else +c itl1=ntortyp +c endif #ifdef MOMENT s1=dip(4,jj,i)*dip(4,kk,k) #endif @@ -9299,7 +9736,7 @@ C energy moment and not to the cluster cumulant. s2=0.5d0*scalar2(b1(1,k),auxvec(1)) call matvec2(AECA(1,1,2),b1(1,l+1),auxvec(1)) s3=0.5d0*scalar2(b1(1,j+1),auxvec(1)) - call transpose2(EE(1,1,itk),auxmat(1,1)) + call transpose2(EE(1,1,k),auxmat(1,1)) call matmat2(auxmat(1,1),AECA(1,1,1),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -9397,25 +9834,25 @@ C C 4/7/01 AL Component s1 was removed, because it pertains to the respective C energy moment and not to the cluster cumulant. cd write (2,*) 'eello_graph4: wturn6',wturn6 - iti=itortyp(itype(i)) - itj=itortyp(itype(j)) - if (j.lt.nres-1) then - itj1=itortyp(itype(j+1)) - else - itj1=ntortyp - endif - itk=itortyp(itype(k)) - if (k.lt.nres-1) then - itk1=itortyp(itype(k+1)) - else - itk1=ntortyp - endif - itl=itortyp(itype(l)) - if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) - else - itl1=ntortyp - endif +c iti=itortyp(itype(i)) +c itj=itortyp(itype(j)) +c if (j.lt.nres-1) then +c itj1=itortyp(itype(j+1)) +c else +c itj1=ntortyp +c endif +c itk=itortyp(itype(k)) +c if (k.lt.nres-1) then +c itk1=itortyp(itype(k+1)) +c else +c itk1=ntortyp +c endif +c itl=itortyp(itype(l)) +c if (l.lt.nres-1) then +c itl1=itortyp(itype(l+1)) +c else +c itl1=ntortyp +c endif cd write (2,*) 'eello6_graph4:','i',i,' j',j,' k',k,' l',l cd write (2,*) 'iti',iti,' itj',itj,' itj1',itj1,' itk',itk, cd & ' itl',itl,' itl1',itl1 @@ -10210,3 +10647,240 @@ C eliptran=elpitran+0.0 ! I am in water enddo return end +C--------------------------------------------------------- +C AFM soubroutine for constant force + subroutine AFMforce(Eafmforce) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.LOCAL' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.NAMES' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.CALC' + include 'COMMON.CONTROL' + include 'COMMON.SPLITELE' + include 'COMMON.SBRIDGE' + real*8 diffafm(3) + dist=0.0d0 + Eafmforce=0.0d0 + do i=1,3 + diffafm(i)=c(i,afmend)-c(i,afmbeg) + dist=dist+diffafm(i)**2 + enddo + dist=dsqrt(dist) + Eafmforce=-forceAFMconst*(dist-distafminit) + do i=1,3 + gradafm(i,afmend-1)=-forceAFMconst*diffafm(i)/dist + gradafm(i,afmbeg-1)=forceAFMconst*diffafm(i)/dist + enddo +C print *,'AFM',Eafmforce + return + end +C--------------------------------------------------------- +C AFM subroutine with pseudoconstant velocity + subroutine AFMvel(Eafmforce) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.LOCAL' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.NAMES' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.CALC' + include 'COMMON.CONTROL' + include 'COMMON.SPLITELE' + include 'COMMON.SBRIDGE' + real*8 diffafm(3) +C Only for check grad COMMENT if not used for checkgrad +C totT=3.0d0 +C-------------------------------------------------------- +C print *,"wchodze" + dist=0.0d0 + Eafmforce=0.0d0 + do i=1,3 + diffafm(i)=c(i,afmend)-c(i,afmbeg) + dist=dist+diffafm(i)**2 + enddo + dist=dsqrt(dist) + Eafmforce=0.5d0*forceAFMconst + & *(distafminit+totTafm*velAFMconst-dist)**2 +C Eafmforce=-forceAFMconst*(dist-distafminit) + do i=1,3 + gradafm(i,afmend-1)=-forceAFMconst* + &(distafminit+totTafm*velAFMconst-dist) + &*diffafm(i)/dist + gradafm(i,afmbeg-1)=forceAFMconst* + &(distafminit+totTafm*velAFMconst-dist) + &*diffafm(i)/dist + enddo +C print *,'AFM',Eafmforce,totTafm*velAFMconst,dist + return + end +C----------------------------------------------------------- +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 +C now the gradient... +C grad_shield is gradient of Calfa for peptide groups + 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 +