iti1=ntortyp+1
endif
cd write (iout,*) '*******i',i,' iti1',iti
-cd write (iout,*) 'b1',b1(:,iti)
-cd write (iout,*) 'b2',b2(:,iti)
+cd write (iout,*) 'b1',b1(:,i-2)
+cd write (iout,*) 'b2',b2(:,i-2)
cd write (iout,*) 'Ug',Ug(:,:,i-2)
c print *,"itilde1 i iti iti1",i,iti,iti1
if (i .gt. iatel_s+2) then
include 'COMMON.IOUNITS'
include 'COMMON.FFIELD'
include 'COMMON.TORCNSTR'
- logical lprn
C Set lprn=.true. for debugging
lprn=.false.
c lprn=.true.
include 'COMMON.IOUNITS'
include 'COMMON.FFIELD'
include 'COMMON.TORCNSTR'
- logical lprn
+ logical lprn,energy_dec
C Set lprn=.true. for debugging
lprn=.false.
c lprn=.true.
+ energy_dec=.false.
etors=0.0D0
do i=iphi_start,iphi_end
if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
itori1=itortyp(itype(i-1))
phii=phi(i)
gloci=0.0D0
+ etors_ii=0.0d0
C Regular cosine and sine terms
do j=1,nterm(itori,itori1,iblock)
v1ij=v1(j,itori,itori1,iblock)
cosphi=dcos(j*phii)
sinphi=dsin(j*phii)
etors=etors+v1ij*cosphi+v2ij*sinphi
+ if (energy_dec) etors_ii=etors_ii+
+ & v1ij*cosphi+v2ij*sinphi
gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
enddo
C Lorentz terms
pom=vl2ij*cosphi+vl3ij*sinphi
pom1=1.0d0/(pom*pom+1.0d0)
etors=etors+vl1ij*pom1
-c if (energy_dec) etors_ii=etors_ii+
-c & vl1ij*pom1
+ if (energy_dec) etors_ii=etors_ii+
+ & vl1ij*pom1
pom=-pom*pom1*pom1
gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom
enddo
C Subtract the constant term
etors=etors-v0(itori,itori1,iblock)
+ if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+ & 'etor',i,etors_ii
if (lprn)
- & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
+ & write (iout,'(2(a3,2x,i3,2x),2i3,f10.2,6f8.3/36x,6f8.3/)')
& restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
+ & rad2deg*phii,
& (v1(j,itori,itori1,1),j=1,6),(v2(j,itori,itori1,1),j=1,6)
gloc(i-3,icg)=gloc(i-3,icg)+wtor*fact*gloci
c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
include 'COMMON.TORCNSTR'
logical lprn
C Set lprn=.true. for debugging
- lprn=.false.
-c lprn=.true.
+c lprn=.false.
+ lprn=.true.
etors_d=0.0D0
do i=iphi_start,iphi_end-1
if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1