X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=inline;f=source%2Funres%2Fsrc_MD%2Fenergy_p_new_barrier.F;h=828e16acdd0e71f6f1df9115457b13e5ec360cd7;hb=664d495e70d14eed4e97f7b8efd2e107dee2fd4e;hp=851ad39228d816dd55e058b3235357e0aa363c6f;hpb=4866688f3356059b9a933ec2be6da24e6784fa9f;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 851ad39..828e16a 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 @@ -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 @@ -8053,9 +8152,9 @@ C C Parallel Antiparallel C C o o -C /l\ /j\ -C / \ / \ -C /| o | | o |\ +C /l\ /j\ +C / \ / \ +C /| o | | o |\ C \ j|/k\| / \ |/k\|l / C \ / \ / \ / \ / C o o o o @@ -8154,18 +8253,18 @@ c---------------------------------------------------------------------------- logical lprn common /kutas/ lprn CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -C -C Parallel Antiparallel -C -C o o -C \ /l\ /j\ / -C \ / \ / \ / -C o| o | | o |o -C \ j|/k\| \ |/k\|l -C \ / \ \ / \ -C o o -C i i -C +C C +C Parallel Antiparallel C +C C +C o o C +C \ /l\ /j\ / C +C \ / \ / \ / C +C o| o | | o |o C +C \ j|/k\| \ |/k\|l C +C \ / \ \ / \ C +C o o C +C i i C +C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cd write (2,*) 'eello6_graph2: i,',i,' j',j,' k',k,' l',l C AL 7/4/01 s1 would occur in the sixth-order moment, @@ -8336,18 +8435,18 @@ c---------------------------------------------------------------------------- double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2) logical swap CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -C -C Parallel Antiparallel -C -C o o -C /l\ / \ /j\ -C / \ / \ / \ -C /| o |o o| o |\ -C j|/k\| / |/k\|l / -C / \ / / \ / -C / o / o -C i i -C +C C +C Parallel Antiparallel C +C C +C o o C +C /l\ / \ /j\ C +C / \ / \ / \ C +C /| o |o o| o |\ C +C j|/k\| / |/k\|l / C +C / \ / / \ / C +C / o / o C +C i i C +C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C 4/7/01 AL Component s1 was removed, because it pertains to the respective @@ -8453,18 +8552,18 @@ c---------------------------------------------------------------------------- & auxvec1(2),auxmat1(2,2) logical swap CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -C -C Parallel Antiparallel -C -C o o -C /l\ / \ /j\ -C / \ / \ / \ -C /| o |o o| o |\ -C \ j|/k\| \ |/k\|l -C \ / \ \ / \ -C o \ o \ -C i i -C +C C +C Parallel Antiparallel C +C C +C o o C +C /l\ / \ /j\ C +C / \ / \ / \ C +C /| o |o o| o |\ C +C \ j|/k\| \ |/k\|l C +C \ / \ \ / \ C +C o \ o \ C +C i i C +C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C 4/7/01 AL Component s1 was removed, because it pertains to the respective