X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Funres%2Fsrc_MD%2Fenergy_p_new_barrier.F;h=097318f301db072b3357603d2a7360d2a10efaac;hb=051632abb45d8a610c7ea1733320282aa1b03587;hp=1a992858b335e7375c346c79b3a4f53a3129ebe4;hpb=478f5fd9e95bd5665ab58ee691f517d7dab379b3;p=unres.git diff --git a/source/unres/src_MD/energy_p_new_barrier.F b/source/unres/src_MD/energy_p_new_barrier.F index 1a99285..097318f 100644 --- a/source/unres/src_MD/energy_p_new_barrier.F +++ b/source/unres/src_MD/energy_p_new_barrier.F @@ -131,6 +131,19 @@ C C Calculate electrostatic (H-bonding) energy of the main chain. C 107 continue + +C BARTEK for dfa test! + if (wdfa_dist.gt.0) call edfad(edfadis) +c print*, 'edfad is finished!', edfadis + if (wdfa_tor.gt.0) call edfat(edfator) +c print*, 'edfat is finished!', edfator + if (wdfa_nei.gt.0) call edfan(edfanei) +c print*, 'edfan is finished!', edfanei + if (wdfa_beta.gt.0) call edfab(edfabet) +c print*, 'edfab is finished!', edfabet +C stop +C BARTEK + c print *,"Processor",myrank," computed USCSC" #ifdef TIMING #ifdef MPI @@ -324,6 +337,10 @@ C energia(21)=esccor energia(22)=evdw_p energia(23)=evdw_m + energia(24)=edfadis + energia(25)=edfator + energia(26)=edfanei + energia(27)=edfabet c print *," Processor",myrank," calls SUM_ENERGY" call sum_energy(energia,.true.) c print *," Processor",myrank," left SUM_ENERGY" @@ -420,6 +437,10 @@ cMS$ATTRIBUTES C :: proc_proc estr=energia(17) Uconst=energia(20) esccor=energia(21) + edfadis=energia(24) + edfator=energia(25) + edfanei=energia(26) + edfabet=energia(27) #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc @@ -427,6 +448,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 + & +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 @@ -434,6 +457,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+Uconst+wsccor*esccor + & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei + & +wdfa_beta*edfabet + #endif energia(0)=etot c detecting NaNQ @@ -540,7 +566,12 @@ c enddo & wcorr5*gradcorr5_long(j,i)+ & wcorr6*gradcorr6_long(j,i)+ & wturn6*gcorr6_turn_long(j,i)+ - & wstrain*ghpbc(j,i) + & wstrain*ghpbc(j,i)+ + & wdfa_dist*gdfad(j,i)+ + & wdfa_tor*gdfat(j,i)+ + & wdfa_nei*gdfan(j,i)+ + & wdfa_beta*gdfab(j,i) + enddo enddo #else @@ -554,7 +585,12 @@ c enddo & wcorr5*gradcorr5_long(j,i)+ & wcorr6*gradcorr6_long(j,i)+ & wturn6*gcorr6_turn_long(j,i)+ - & wstrain*ghpbc(j,i) + & wstrain*ghpbc(j,i)+ + & wdfa_dist*gdfad(j,i)+ + & wdfa_tor*gdfat(j,i)+ + & wdfa_nei*gdfan(j,i)+ + & wdfa_beta*gdfab(j,i) + enddo enddo #endif @@ -570,7 +606,13 @@ c enddo & wcorr5*gradcorr5_long(j,i)+ & wcorr6*gradcorr6_long(j,i)+ & wturn6*gcorr6_turn_long(j,i)+ - & wstrain*ghpbc(j,i) + & wstrain*ghpbc(j,i)+ + & wdfa_dist*gdfad(j,i)+ + & wdfa_tor*gdfat(j,i)+ + & wdfa_nei*gdfan(j,i)+ + & wdfa_beta*gdfab(j,i) + + enddo enddo #endif @@ -756,6 +798,7 @@ c enddo & +wturn3*gel_loc_turn3(i) & +wturn6*gel_loc_turn6(i) & +wel_loc*gel_loc_loc(i) + & +wsccor*gsccor_loc(i) enddo #ifdef DEBUG write (iout,*) "gloc after adding corr" @@ -1046,6 +1089,12 @@ C------------------------------------------------------------------------ estr=energia(17) Uconst=energia(20) esccor=energia(21) +C Bartek + edfadis = energia(24) + edfator = energia(25) + edfanei = energia(26) + edfabet = energia(27) + #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp, & estr,wbond,ebe,wang, @@ -1054,7 +1103,7 @@ C------------------------------------------------------------------------ & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor, & edihcnstr,ebr*nss, - & Uconst,etot + & Uconst,edfadis,edfator,edfanei,edfabet,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -1065,7 +1114,7 @@ C------------------------------------------------------------------------ & '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, + & 'EHPB= ',1pE16.6,' WEIGHT=',1pD16.6, & ' (SS bridges & dist. cnstr.)'/ & 'ECORR4=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ & 'ECORR5=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ @@ -1078,6 +1127,10 @@ C------------------------------------------------------------------------ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST= ',1pE16.6,' (Constraint energy)'/ + & 'EDFAD= ',1pE16.6,' (DFA distance energy)'/ + & 'EDFAT= ',1pE16.6,' (DFA torsion energy)'/ + & 'EDFAN= ',1pE16.6,' (DFA NCa energy)'/ + & 'EDFAB= ',1pE16.6,' (DFA Beta energy)'/ & 'ETOT= ',1pE16.6,' (total)') #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec, @@ -1086,7 +1139,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,etot + & ebr*nss, + & Uconst,edfadis,edfator,edfanei,edfabet,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -1109,6 +1163,10 @@ C------------------------------------------------------------------------ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST=',1pE16.6,' (Constraint energy)'/ + & 'EDFAD= ',1pE16.6,' (DFA distance energy)'/ + & 'EDFAT= ',1pE16.6,' (DFA torsion energy)'/ + & 'EDFAN= ',1pE16.6,' (DFA NCa energy)'/ + & 'EDFAB= ',1pE16.6,' (DFA Beta energy)'/ & 'ETOT= ',1pE16.6,' (total)') #endif return @@ -1447,7 +1505,7 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) chi1=chi(itypi,itypj) @@ -4266,49 +4324,90 @@ C iii and jjj point to the residues for which the distance is assigned. iii=ii jjj=jj endif -cd write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj +c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj, +c & dhpb(i),dhpb1(i),forcon(i) C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij cd write (iout,*) "eij",eij + else if (ii.gt.nres .and. jj.gt.nres) then +c Restraints from contact prediction + dd=dist(ii,jj) + 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 + 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) - rdis=dd-dhpb(i) + dd=dist(ii,jj) + 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) + waga=forcon(i) C Calculate the contribution to energy. - ehpb=ehpb+waga*rdis*rdis + ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "alpha reg",dd,waga*rdis*rdis C C Evaluate gradient. C - fac=waga*rdis/dd + fac=waga*rdis/dd + endif cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, cd & ' waga=',waga,' fac=',fac - do j=1,3 - ggg(j)=fac*(c(j,jj)-c(j,ii)) - enddo + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3) C If this is a SC-SC distance, we need to calculate the contributions to the C Cartesian gradient in the SC vectors (ghpbx). - if (iii.lt.ii) then + if (iii.lt.ii) then do j=1,3 ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) enddo - endif + endif cgrad do j=iii,jjj-1 cgrad do k=1,3 cgrad ghpbc(k,j)=ghpbc(k,j)+ggg(k) cgrad enddo cgrad enddo - do k=1,3 - ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) - ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) - enddo + do k=1,3 + ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) + ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) + enddo endif enddo ehpb=0.5D0*ehpb @@ -5254,7 +5353,7 @@ c do j = 1,3 xx = xx + x_prime(j)*dc_norm(j,i+nres) yy = yy + y_prime(j)*dc_norm(j,i+nres) - zz = zz + z_prime(j)*dc_norm(j,i+nres) + zz = zz + dsign(1.0,itype(i))*z_prime(j)*dc_norm(j,i+nres) enddo xxtab(i)=xx @@ -5264,7 +5363,7 @@ C C Compute the energy of the ith side cbain C c write (2,*) "xx",xx," yy",yy," zz",zz - it=itype(i) + it=iabs(itype(i)) do j = 1,65 x(j) = sc_parmin(j,it) enddo @@ -5272,7 +5371,7 @@ c write (2,*) "xx",xx," yy",yy," zz",zz Cc diagnostics - remove later xx1 = dcos(alph(2)) yy1 = dsin(alph(2))*dcos(omeg(2)) - zz1 = -dsin(alph(2))*dsin(omeg(2)) + zz1 = -dsign(1.0,itype(i))*dsin(alph(2))*dsin(omeg(2)) write(2,'(3f8.1,3f9.3,1x,3f9.3)') & alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz, & xx1,yy1,zz1 @@ -5898,6 +5997,7 @@ c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor esccor=0.0D0 do i=itau_start,itau_end esccor_ii=0.0D0 + if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle isccori=isccortyp(itype(i-2)) isccori1=isccortyp(itype(i-1)) phii=phi(i) @@ -5916,14 +6016,14 @@ c 2 = Ca...Ca...Ca...SC c 3 = SC...Ca...Ca...SCi gloci=0.0D0 if (((intertyp.eq.3).and.((itype(i-2).eq.10).or. - & (itype(i-1).eq.10).or.(itype(i-2).eq.21).or. - & (itype(i-1).eq.21))) + & (itype(i-1).eq.10).or.(itype(i-2).eq.ntyp1).or. + & (itype(i-1).eq.ntyp1))) & .or. ((intertyp.eq.1).and.((itype(i-2).eq.10) - & .or.(itype(i-2).eq.21))) + & .or.(itype(i-2).eq.ntyp1))) & .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or. - & (itype(i-1).eq.21)))) cycle - if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.21)) cycle - if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.21)) + & (itype(i-1).eq.ntyp1)))) cycle + if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.ntyp1)) cycle + if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.ntyp1)) & cycle do j=1,nterm_sccor(isccori,isccori1) v1ij=v1sccor(j,intertyp,isccori,isccori1)