X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fenergy_p_new_barrier.F;h=7cf67dba1c3afaa0bc6a74f506e8e8fd5c272915;hb=a4e5a256e77beab0f13a18e273b46707d587060a;hp=9edadf8492003c7338bead2cbda03c0ad18892c6;hpb=52187c9365594597457799445140ea830f1141a3;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 9edadf8..7cf67db 100644 --- a/source/unres/src_MD/energy_p_new_barrier.F +++ b/source/unres/src_MD/energy_p_new_barrier.F @@ -27,6 +27,7 @@ cMS$ATTRIBUTES C :: proc_proc #ifdef MPI c print*,"ETOTAL Processor",fg_rank," absolute rank",myrank, c & " nfgtasks",nfgtasks + call flush(iout) if (nfgtasks.gt.1) then #ifdef MPI time00=MPI_Wtime() @@ -131,6 +132,36 @@ C C Calculate electrostatic (H-bonding) energy of the main chain. C 107 continue +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 +c print*, 'edfab is finished!', edfabet +cmc +cmc Sep-06: egb takes care of dynamic ss bonds too +cmc +c if (dyn_ss) call dyn_set_nss + c print *,"Processor",myrank," computed USCSC" #ifdef TIMING #ifdef MPI @@ -223,6 +254,15 @@ cd print *,'nterm=',nterm etors=0 edihcnstr=0 endif + + if (constr_homology.ge.1) then + call e_modeller(ehomology_constr) + else + ehomology_constr=0.0d0 + endif + + +c write(iout,*) ehomology_constr c print *,"Processor",myrank," computed Utor" C C 6/23/01 Calculate double-torsional energy @@ -324,8 +364,14 @@ C energia(21)=esccor energia(22)=evdw_p energia(23)=evdw_m + energia(24)=ehomology_constr + energia(25)=edfadis + energia(26)=edfator + energia(27)=edfanei + energia(28)=edfabet c print *," Processor",myrank," calls SUM_ENERGY" call sum_energy(energia,.true.) + if (dyn_ss) call dyn_set_nss c print *," Processor",myrank," left SUM_ENERGY" #ifdef TIMING #ifdef MPI @@ -420,20 +466,29 @@ cMS$ATTRIBUTES C :: proc_proc estr=energia(17) Uconst=energia(20) esccor=energia(21) + ehomology_constr=energia(24) + edfadis=energia(25) + edfator=energia(26) + edfanei=energia(27) + edfabet=energia(28) #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 + & +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 + & +wbond*estr+Uconst+wsccor*esccor+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 - & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 + & +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 + & +wbond*estr+Uconst+wsccor*esccor+ehomology_constr + & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei + & +wdfa_beta*edfabet #endif energia(0)=etot c detecting NaNQ @@ -540,7 +595,11 @@ 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 +613,11 @@ 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 +633,11 @@ 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 @@ -774,7 +841,6 @@ c enddo do i=1,4*nres glocbuf(i)=gloc(i,icg) enddo -#define DEBUG #ifdef DEBUG write (iout,*) "gloc_sc before reduce" do i=1,nres @@ -783,7 +849,6 @@ c enddo enddo enddo #endif -#undef DEBUG do i=1,nres do j=1,3 gloc_scbuf(j,i)=gloc_sc(j,i,icg) @@ -802,7 +867,6 @@ c enddo call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres, & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) time_reduce=time_reduce+MPI_Wtime()-time00 -#define DEBUG #ifdef DEBUG write (iout,*) "gloc_sc after reduce" do i=1,nres @@ -811,7 +875,6 @@ c enddo enddo enddo #endif -#undef DEBUG #ifdef DEBUG write (iout,*) "gloc after reduce" do i=1,4*nres @@ -1046,6 +1109,13 @@ C------------------------------------------------------------------------ estr=energia(17) Uconst=energia(20) esccor=energia(21) + ehomology_constr=energia(24) +C Bartek + edfadis = energia(25) + edfator = energia(26) + edfanei = energia(27) + edfabet = energia(28) + #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp, & estr,wbond,ebe,wang, @@ -1053,31 +1123,37 @@ C------------------------------------------------------------------------ & 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,etot + & edihcnstr,ehomology_constr, ebr*nss, + & Uconst,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)'/ - & 'EHPB= ',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)'/ + & 'EHPB= ',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)'/ + & 'ECORR4=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ + & 'ECORR5=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ + & 'ECORR6=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ + & '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)'/ + & 'ETURN6=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 6th order)'/ + & 'ESCCOR=',1pE16.6,' WEIGHT=',1pE16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST= ',1pE16.6,' (Constraint 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, @@ -1086,7 +1162,9 @@ 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 + & ehomology_constr,ebr*nss,Uconst,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)'/ @@ -1107,8 +1185,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)'/ + & 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST=',1pE16.6,' (Constraint energy)'/ + & 'EDFAD= ',1pE16.6,' WEIGHT=',1pD16.6,' (DFA distance energy)'/ + & 'EDFAT= ',1pE16.6,' WEIGHT=',1pD16.6,' (DFA torsion energy)'/ + & 'EDFAN= ',1pE16.6,' WEIGHT=',1pD16.6,' (DFA NCa energy)'/ + & 'EDFAB= ',1pE16.6,' WEIGHT=',1pD16.6,' (DFA Beta energy)'/ & 'ETOT= ',1pE16.6,' (total)') #endif return @@ -1556,6 +1639,7 @@ C include 'COMMON.IOUNITS' include 'COMMON.CALC' include 'COMMON.CONTROL' + include 'COMMON.SBRIDGE' logical lprn evdw=0.0D0 ccccc energy_dec=.false. @@ -1584,6 +1668,12 @@ C Calculate SC interaction energy. 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 + call dyn_ssbond_ene(i,j,evdwij) + evdw=evdw+evdwij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') + & 'evdw',i,j,evdwij,' ss' + ELSE ind=ind+1 itypj=itype(j) c dscj_inv=dsc_inv(itypj) @@ -1693,6 +1783,7 @@ C Calculate angular part of the gradient. #else call sc_grad #endif + ENDIF ! dyn_ss enddo ! j enddo ! iint enddo ! i @@ -4270,9 +4361,15 @@ 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 +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 call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij + endif cd write (iout,*) "eij",eij else if (ii.gt.nres .and. jj.gt.nres) then c Restraints from contact prediction @@ -4411,7 +4508,7 @@ c dscj_inv=dsc_inv(itypj) deltat12=om2-om1+2.0d0 cosphi=om12-om1*om2 eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2) - & +akct*deltad*deltat12 + & +akct*deltad*deltat12+ebr & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi c write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth, c & " akct",akct," deltad",deltad," deltat",deltat1,deltat2, @@ -4621,7 +4718,7 @@ C Derivatives of the "mean" values in gamma1 and gamma2. & 'ebend',i,ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1 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)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg) enddo C Ufff.... We've done all this!!! return @@ -4912,13 +5009,15 @@ C enddo enddo 10 continue - if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') - & i,theta(i)*rad2deg,phii*rad2deg, +c lprn1=.true. + if (lprn1) write (iout,'(a4,i2,3f8.1,9h ethetai ,f10.5)') + & 'ebe', i,theta(i)*rad2deg,phii*rad2deg, & phii1*rad2deg,ethetai +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)+wang*dethetai enddo return end @@ -5742,6 +5841,15 @@ C Proline-Proline pair is a special case... return end c------------------------------------------------------------------------------ +c LICZENIE WIEZOW Z ROWNANIA ENERGII MODELLERA + subroutine e_modeller(ehomology_constr) + ehomology_constr=0.0 + write (iout,*) "!!!!!UWAGA, JESTEM W DZIWNEJ PETLI, TEST!!!!!" + return + end +C !!!!!!!! NIE CZYTANE !!!!!!!!!!! + +c------------------------------------------------------------------------------ subroutine etor_d(etors_d) etors_d=0.0d0 return @@ -5842,6 +5950,187 @@ cd write (iout,*) 'edihcnstr',edihcnstr return end c---------------------------------------------------------------------------- +c MODELLER restraint function + subroutine e_modeller(ehomology_constr) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + + integer nnn, i, j, k, ki, irec, l + integer katy, odleglosci, test7 + real*8 odleg, odleg2, odleg3, kat, kat2, kat3, gdih(max_template) + real*8 distance(max_template),distancek(max_template), + & min_odl,godl(max_template),dih_diff(max_template) + + include 'COMMON.SBRIDGE' + include 'COMMON.CHAIN' + include 'COMMON.GEO' + include 'COMMON.DERIV' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.VAR' + include 'COMMON.IOUNITS' + include 'COMMON.MD' + include 'COMMON.CONTROL' + + + do i=1,19 + 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 + do ii = link_start_homo,link_end_homo + i = ires_homo(ii) + j = jres_homo(ii) + dij=dist(i,j) + do k=1,constr_homology + distance(k)=odl(k,ii)-dij + distancek(k)= + & 0.5d0*waga_dist(iset)*distance(k)**2*sigma_odl(k,ii) + enddo + + min_odl=minval(distancek) +#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 + odleg2=0.0d0 + 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)) + godl(k)=dexp(-distancek(k)+min_odl) + odleg2=odleg2+godl(k) + +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 +#ifdef DEBUG + write (iout,*) "godl",(godl(k),k=1,constr_homology) + write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 +#endif + odleg=odleg-dLOG(odleg2/constr_homology)+min_odl +c Gradient + sum_godl=odleg2 + sum_sgodl=0.0 + do k=1,constr_homology +c godl=dexp(((-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2)) +c & *waga_dist(iset))+min_odl + sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist(iset) + 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 + + 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) + + enddo +ccc write(iout,778)"TEST: odleg2=", odleg2, "DLOG(odleg2)=", +ccc & dLOG(odleg2),"-odleg=", -odleg + + enddo ! ii +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 + do i=idihconstr_start_homo,idihconstr_end_homo + kat2=0.0d0 +c betai=beta(i,i+1,i+2,i+3) + betai = phi(i+3) + do k=1,constr_homology + dih_diff(k)=pinorm(dih(k,i)-betai) +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) + + kat3=-0.5d0*waga_angle(iset)*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 +#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) + +ccc write(iout,778)"TEST: kat2=", kat2, "DLOG(kat2)=", +ccc & dLOG(kat2), "-kat=", -kat + +c ---------------------------------------------------------------------- +c Gradient +c ---------------------------------------------------------------------- + + sum_gdih=kat2 + sum_sgdih=0.0 + do k=1,constr_homology + sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle(iset) + sum_sgdih=sum_sgdih+sgdih + enddo + grad_dih3=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,icg)=gloc(i,icg)+grad_dih3 +ccc write(iout,747) "GRAD_KAT_2", i, nphi, icg, grad_dih3, +ccc & gloc(nphi+i-3,icg) + + enddo + + +c Total energy from homology restraints +#ifdef DEBUG + write (iout,*) "odleg",odleg," kat",kat +#endif + ehomology_constr=odleg+kat + return + + 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------------------------------------------------------------------------------ subroutine etor_d(etors_d) C 6/23/01 Compute double torsional energy implicit real*8 (a-h,o-z) @@ -8191,7 +8480,7 @@ c---------------------------------------------------------------------------- include 'COMMON.GEO' logical swap double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2), - & auxvec1(2),auxvec2(1),auxmat1(2,2) + & auxvec1(2),auxvec2(2),auxmat1(2,2) logical lprn common /kutas/ lprn CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC