X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fenergy_p_new_barrier.F;h=cbf930c6d4bcb460f0174fe5ccedb9ad1ba985da;hb=f37cfff89ef75a184b0c7bbbac2ecea2efc1d157;hp=535e8b2ae59729c9862a7fa7d3333b2c46b52682;hpb=c45f3ac598572d95b62213f9ed78f4990a8f7735;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 535e8b2..cbf930c 100644 --- a/source/unres/src_MD/energy_p_new_barrier.F +++ b/source/unres/src_MD/energy_p_new_barrier.F @@ -99,6 +99,12 @@ c if (modecalc.eq.12.or.modecalc.eq.14) then c call int_from_cart1(.false.) c endif #endif +#ifndef DFA + edfadis=0.0d0 + edfator=0.0d0 + edfanei=0.0d0 + edfabet=0.0d0 +#endif #ifdef TIMING #ifdef MPI time00=MPI_Wtime() @@ -132,6 +138,7 @@ 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) @@ -156,6 +163,7 @@ c print*, 'edfan is finished!', edfanei else edfabet=0 endif +#endif c print*, 'edfab is finished!', edfabet cmc cmc Sep-06: egb takes care of dynamic ss bonds too @@ -255,11 +263,11 @@ cd print *,'nterm=',nterm edihcnstr=0 endif - if (constr_homology.ge.1) then + if (constr_homology.ge.1.and.waga_homology(iset).ne.0d0) then call e_modeller(ehomology_constr) - print *,'iset=',iset,'me=',me,ehomology_constr, - & 'Processor',fg_rank,' CG group',kolor, - & ' absolute rank',MyRank +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 @@ -543,6 +551,7 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.TIME1' include 'COMMON.MAXGRAD' include 'COMMON.SCCOR' + include 'COMMON.MD' #ifdef TIMING #ifdef MPI time01=MPI_Wtime() @@ -813,6 +822,14 @@ c enddo #endif enddo enddo + if (constr_homology.gt.0.and.waga_homology(iset).ne.0d0) 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,*) "gloc before adding corr" do i=1,4*nres @@ -1764,9 +1781,10 @@ 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) then + write (iout,'(a6,2i5,0pf7.3)') 'evdw',i,j,evdwij + call flush(iout) + endif C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 fac=-expon*(e1+evdwij)*rij_shift @@ -3877,6 +3895,7 @@ C Derivatives in gamma(i+1) gel_loc_turn3(i+1)=gel_loc_turn3(i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) C Cartesian derivatives +!DIR$ UNROLL(0) do l=1,3 c ghalf1=0.5d0*agg(l,1) c ghalf2=0.5d0*agg(l,2) @@ -4343,6 +4362,7 @@ C include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' + include 'COMMON.CONTROL' dimension ggg(3) ehpb=0.0D0 cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr @@ -4378,6 +4398,16 @@ 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 (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,f15.6,2f8.3)') + & "edisl",ii,jj, + & fordepth(i)**4.0d0*rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)), + & 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 @@ -4395,7 +4425,8 @@ C C Evaluate gradient. C fac=waga*rdis/dd - endif + endif + endif do j=1,3 ggg(j)=fac*(c(j,jj)-c(j,ii)) enddo @@ -4411,6 +4442,18 @@ C 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,f15.6,2f8.3)') + 7 "edisl",ii,jj, + & fordepth(i)**4.0d0*rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)), + & fordepth(i),dd +c if (energy_dec) +c & write (iout,*) fac + 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 @@ -4428,6 +4471,7 @@ C Evaluate gradient. C fac=waga*rdis/dd endif + endif cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, cd & ' waga=',waga,' fac=',fac do j=1,3 @@ -4453,7 +4497,10 @@ cgrad enddo enddo endif enddo - ehpb=0.5D0*ehpb + if (constr_dist.ne.11) ehpb=0.5D0*ehpb +c do i=1,nres +c write (iout,*) "ghpbc",i,(ghpbc(j,i),j=1,3) +c enddo return end C-------------------------------------------------------------------------- @@ -4567,6 +4614,8 @@ c do i=ibondp_start,ibondp_end diff = vbld(i)-vbldp0 c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff + 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 gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i) @@ -4585,6 +4634,12 @@ c diff=vbld(i+nres)-vbldsc0(1,iti) c write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff, c & AKSC(1,iti),AKSC(1,iti)*diff*diff + if (energy_dec) then + write (iout,*) + & "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff, + & AKSC(1,iti),AKSC(1,iti)*diff*diff + call flush(iout) + endif estr=estr+0.5d0*AKSC(1,iti)*diff*diff do j=1,3 gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres) @@ -4862,6 +4917,8 @@ C & sinph1ph2(maxdouble,maxdouble) logical lprn /.false./, lprn1 /.false./ etheta=0.0D0 +c write (iout,*) "EBEND ithet_start",ithet_start, +c & " ithet_end",ithet_end do i=ithet_start,ithet_end if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. &(itype(i).eq.ntyp1)) cycle @@ -4875,7 +4932,7 @@ C sinkt(k)=dsin(k*theti2) enddo C if (i.gt.3) then - if (i.gt.3 .and. itype(i-3).ne.ntyp1) then + if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -5022,6 +5079,8 @@ c lprn1=.true. & phii1*rad2deg,ethetai c lprn1=.false. etheta=etheta+ethetai + if (energy_dec) write (iout,'(a6,i5,0pf7.3)') + & 'ebend',i,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)=gloc(nphi+i-2,icg)+wang*dethetai @@ -5851,7 +5910,7 @@ C Proline-Proline pair is a special case... c------------------------------------------------------------------------------ c LICZENIE WIEZOW Z ROWNANIA ENERGII MODELLERA subroutine e_modeller(ehomology_constr) - ehomology_constr=0.0 + ehomology_constr=0.0d0 write (iout,*) "!!!!!UWAGA, JESTEM W DZIWNEJ PETLI, TEST!!!!!" return end @@ -5999,7 +6058,7 @@ c c - do i=1,19 + do i=1,max_template distancek(i)=9999999.9 enddo @@ -6018,16 +6077,45 @@ c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d 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 write (iout,*) "distance: ii",ii," nexl",nexl - min_odl=minval(distancek) + +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 @@ -6035,13 +6123,32 @@ c write (iout,* )"min_odl",min_odl 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=", @@ -6055,16 +6162,42 @@ c write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps write (iout,*) "godl",(godl(k),k=1,constr_homology) ! exponents write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps #endif - odleg=odleg-dLOG(odleg2/constr_homology)+min_odl + 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 - sum_godl=odleg2 +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 @@ -6072,12 +6205,22 @@ 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 (homol_nset.gt.1)then - grad_odl3=waga_dist1(iset)*sum_sgodl/(sum_godl*dij) - else - grad_odl3=waga_dist*sum_sgodl/(sum_godl*dij) - endif + 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) @@ -6131,10 +6274,10 @@ c write (iout,*) idihconstr_start_homo,idihconstr_end_homo 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+3) + betai = phi(i) c write (iout,*) "betai =",betai + kat2=0.0d0 do k=1,constr_homology dih_diff(k)=pinorm(dih(k,i)-betai) c write (iout,*) "dih_diff(",k,") =",dih_diff(k) @@ -6142,13 +6285,16 @@ 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(*,*)"" +c write(iout,*) "i",i," k",k," sigma",sigma_dih(k,i), +c & " kat2=", kat2, "gdih=",gdih(k) enddo c write (iout,*) "gdih",(gdih(k),k=1,constr_homology) ! exps c write (iout,*) "i",i," betai",betai," kat2",kat2 ! sum of exps @@ -6170,21 +6316,21 @@ 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 - if (homol_nset.gt.1)then - grad_dih3=waga_angle1(iset)*sum_sgdih/sum_gdih - else - grad_dih3=waga_angle*sum_sgdih/sum_gdih - endif + 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,icg)=gloc(i,icg)+grad_dih3 + 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 @@ -6253,7 +6399,9 @@ 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+dexp(utheta_i) ! Sum of Gaussians (pk) + gutheta_i=gutheta_i+gtheta(k) ! Sum of Gaussians (pk) +c write (iout,*) "i",i," k",k," sigma_theta",sigma_theta(k,i), +c & " gtheta",gtheta(k) 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 @@ -6272,11 +6420,12 @@ 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 grad_theta3=sum_sgtheta/sum_gtheta 1/*theta(i)? s. line below -c grad_theta3=sum_sgtheta/sum_gtheta -c c Final value of gradient using same var as in Econstr_back - dutheta(i-2)=sum_sgtheta/sum_gtheta*waga_theta + 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 @@ -6324,7 +6473,9 @@ c 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 - guscdiff(i)=guscdiff(i)+dexp(usc_diff_i) !Sum of Gaussians (pk) +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 @@ -6365,7 +6516,7 @@ c sum_sguscdiff=sum_sguscdiff+sum_guscdiff c c c New implementation - sum_guscdiff = waga_d*sum_guscdiff + 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+ @@ -6453,16 +6604,31 @@ c c Addition of energy of theta angle and SC local geom over constr_homologs ref strs c c ehomology_constr=odleg+kat - if (homol_nset.gt.1)then - ehomology_constr=waga_dist1(iset)*odleg+waga_angle1(iset)*kat+waga_theta*Eval - & +waga_d*Erot +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 - ehomology_constr=waga_dist*odleg+waga_angle*kat+waga_theta*Eval - & +waga_d*Erot +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 -c write (iout,*) "odleg",odleg," kat",kat," Uconst_back",Uconst_back -c write (iout,*) "ehomology_constr",ehomology_constr -c ehomology_constr=odleg+kat+Uconst_back +#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 @@ -6491,12 +6657,14 @@ C 6/23/01 Compute double torsional energy include 'COMMON.IOUNITS' include 'COMMON.FFIELD' include 'COMMON.TORCNSTR' + include 'COMMON.CONTROL' logical lprn C Set lprn=.true. for debugging lprn=.false. c lprn=.true. etors_d=0.0D0 do i=iphid_start,iphid_end + etors_d_ii=0.0D0 itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) itori2=itortyp(itype(i)) @@ -6515,6 +6683,8 @@ c lprn=.true. sinphi2=dsin(j*phii1) etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+ & v2cij*cosphi2+v2sij*sinphi2 + if (energy_dec) etors_d_ii=etors_d_ii+ + & v1cij*cosphi1+v1sij*sinphi1+v2cij*cosphi2+v2sij*sinphi2 gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1) gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2) enddo @@ -6530,12 +6700,17 @@ c lprn=.true. sinphi1m2=dsin(l*phii-(k-l)*phii1) etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+ & v1sdij*sinphi1p2+v2sdij*sinphi1m2 + if (energy_dec) etors_d_ii=etors_d_ii+ + & v1cdij*cosphi1p2+v2cdij*cosphi1m2+ + & v1sdij*sinphi1p2+v2sdij*sinphi1m2 gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2 & -v1cdij*sinphi1p2-v2cdij*sinphi1m2) gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2 & -v1cdij*sinphi1p2+v2cdij*sinphi1m2) enddo enddo + if (energy_dec) write (iout,'(a6,i5,0pf7.3)') + & 'etor_d',i,etors_d_ii gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1 gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2 c write (iout,*) "gloci", gloc(i-3,icg)