C Calculate electrostatic (H-bonding) energy of the main chain.
C
107 continue
+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
time01=MPI_Wtime()
energia(17)=estr
energia(20)=Uconst+Uconst_back
energia(21)=esccor
+c Here are the energies showed per procesor if the are more processors
+c per molecule then we sum it up in sum_energy subroutine
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
time_sumene=time_sumene+MPI_Wtime()-time00
#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
#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
do i=1,4*nres
glocbuf(i)=gloc(i,icg)
enddo
-#define DEBUG
+c#define DEBUG
#ifdef DEBUG
write (iout,*) "gloc_sc before reduce"
do i=1,nres
enddo
enddo
#endif
-#undef DEBUG
+c#undef DEBUG
do i=1,nres
do j=1,3
gloc_scbuf(j,i)=gloc_sc(j,i,icg)
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
+c#define DEBUG
#ifdef DEBUG
write (iout,*) "gloc_sc after reduce"
do i=1,nres
enddo
enddo
#endif
-#undef DEBUG
+c#undef DEBUG
#ifdef DEBUG
write (iout,*) "gloc after reduce"
do i=1,4*nres
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
include 'COMMON.CONTROL'
+ include 'COMMON.SBRIDGE'
logical lprn
+
+c write(iout,*) "Jestem w egb(evdw)"
+
evdw=0.0D0
ccccc energy_dec=.false.
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
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
+
+c write(iout,*) "PRZED ZWYKLE", evdwij
+ call dyn_ssbond_ene(i,j,evdwij)
+c write(iout,*) "PO ZWYKLE", evdwij
+
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+ & 'evdw',i,j,evdwij,' ss'
+C triple bond artifac removal
+ do k=j+1,iend(i,iint)
+C search over all next residues
+ if (dyn_ss_mask(k)) then
+C check if they are cysteins
+C write(iout,*) 'k=',k
+
+c write(iout,*) "PRZED TRI", evdwij
+ evdwij_przed_tri=evdwij
+ call triple_ssbond_ene(i,j,k,evdwij)
+c if(evdwij_przed_tri.ne.evdwij) then
+c write (iout,*) "TRI:", evdwij, evdwij_przed_tri
+c endif
+
+c write(iout,*) "PO TRI", evdwij
+C call the energy function that removes the artifical triple disulfide
+C bond the soubroutine is located in ssMD.F
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+ & 'evdw',i,j,evdwij,'tss'
+ endif!dyn_ss_mask(k)
+ enddo! k
+ ELSE
ind=ind+1
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
gg(3)=zj*fac
C Calculate angular part of the gradient.
call sc_grad
+ ENDIF ! dyn_ss
enddo ! j
enddo ! iint
enddo ! i
enddo
c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
if (i.gt. nnt+1 .and. i.lt.nct+1) then
- iti1 = itortyp(itype(i-1))
+ if (itype(i-1).le.ntyp) then
+ iti1 = itortyp(itype(i-1))
+ else
+ iti1=ntortyp+1
+ endif
else
iti1=ntortyp+1
endif
cd & xmedi,ymedi,zmedi,xj,yj,zj
if (energy_dec) then
- write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
+ write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)')
+ &'evdw1',i,j,evdwij
+ &,iteli,itelj,aaa,evdw1
write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
endif
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
& 'eelloc',i,j,eel_loc_ij
+c write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
eel_loc=eel_loc+eel_loc_ij
C Partial derivatives in virtual-bond dihedral angles gamma
endif
evdwij=e1+e2
evdw2=evdw2+evdwij
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
- & 'evdw2',i,j,evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)')
+ & 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),
+ & bad(itypj,iteli)
C
C Calculate contributions to the gradient in the virtual-bond and SC vectors.
C
include 'COMMON.VAR'
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
+ include 'COMMON.CONTROL'
dimension ggg(3)
ehpb=0.0D0
+ do i=1,3
+ ggg(i)=0.0d0
+ enddo
+C write (iout,*) ,"link_end",link_end,constr_dist
cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
cd write(iout,*)'link_start=',link_start,' link_end=',link_end
if (link_end.eq.0) return
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. iabs(itype(iii)).eq.1 .and.
+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. iabs(itype(iii)).eq.1 .and.
& iabs(itype(jjj)).eq.1) then
call ssbond_ene(iii,jjj,eij)
ehpb=ehpb+2*eij
+ endif
cd write (iout,*) "eij",eij
+cd & ' waga=',waga,' fac=',fac
+ 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,3f8.3)') "edisl",ii,jj,
+ & ehpb,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
+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
+ 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 (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,3f8.3)') "edisl",ii,jj,
+ & ehpb,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
+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
-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
+ fac=waga*rdis/dd
+ endif
+ endif
+ 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
+ if (constr_dist.ne.11) ehpb=0.5D0*ehpb
return
end
C--------------------------------------------------------------------------
cosphi=om12-om1*om2
eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2)
& +akct*deltad*deltat12
- & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
+ & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi+ebr
c write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth,
c & " akct",akct," deltad",deltad," deltat",deltat1,deltat2,
c & " deltat12",deltat12," eij",eij
logical lprn /.false./, lprn1 /.false./
etheta=0.0D0
do i=ithet_start,ithet_end
- if (itype(i-1).eq.ntyp1) cycle
+ if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+ &(itype(i).eq.ntyp1)) cycle
+C print *,i,theta(i)
if (iabs(itype(i+1)).eq.20) iblock=2
if (iabs(itype(i+1)).ne.20) iblock=1
dethetai=0.0d0
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
- if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+C print *,ethetai
+
+ if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
enddo
else
phii=0.0d0
- ityp1=nthetyp+1
do k=1,nsingle
+ ityp1=ithetyp((itype(i-2)))
cosph1(k)=0.0d0
sinph1(k)=0.0d0
enddo
endif
- if (i.lt.nres .and. itype(i).ne.ntyp1) then
+ if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
if (phii1.ne.phii1) phii1=150.0
enddo
else
phii1=0.0d0
- ityp3=nthetyp+1
+ ityp3=ithetyp((itype(i)))
do k=1,nsingle
cosph2(k)=0.0d0
sinph2(k)=0.0d0
sinph1ph2(k,l)=scl-csl
enddo
enddo
-c if (lprn) then
+ if (lprn) then
write (iout,*) "i",i," ityp1",ityp1," ityp2",ityp2,
& " ityp3",ityp3," theti2",theti2," phii",phii," phii1",phii1
-c write (iout,*) "coskt and sinkt"
-c do k=1,nntheterm
-c write (iout,*) k,coskt(k),sinkt(k)
-c enddo
-c endif
+ write (iout,*) "coskt and sinkt"
+ do k=1,nntheterm
+ write (iout,*) k,coskt(k),sinkt(k)
+ enddo
+ endif
do k=1,ntheterm
ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3,iblock)*sinkt(k)
dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3,iblock)
enddo
write(iout,*) "ethetai",ethetai
endif
+C print *,ethetai
do m=1,ntheterm2
do k=1,nsingle
aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)
& ccthet(k,m,ityp1,ityp2,ityp3,iblock)," ddthet",
& ddthet(k,m,ityp1,ityp2,ityp3,iblock)," eethet",
& eethet(k,m,ityp1,ityp2,ityp3,iblock)," ethetai",ethetai
+C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k)
enddo
enddo
+C print *,"cosph1", (cosph1(k), k=1,nsingle)
+C print *,"cosph2", (cosph2(k), k=1,nsingle)
+C print *,"sinph1", (sinph1(k), k=1,nsingle)
+C print *,"sinph2", (sinph2(k), k=1,nsingle)
if (lprn)
& write(iout,*) "ethetai",ethetai
+C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k)
do m=1,ntheterm3
do k=2,ndouble
do l=1,k-1
enddo
10 continue
c lprn1=.true.
-c if (lprn1)
- write (iout,'(i2,3f8.1,9h ethetai ,f10.5)')
+C print *,ethetai
+ if (lprn1)
+ & write (iout,'(i2,3f8.1,9h ethetai ,f10.5)')
& 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
cosfac=dsqrt(cosfac2)
sinfac2=0.5d0/(1.0d0-costtab(i+1))
sinfac=dsqrt(sinfac2)
- it=itype(i)
+ it=iabs(itype(i))
if (it.eq.10) goto 1
c
C Compute the axes of tghe local cartesian coordinates system; store in
y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
enddo
do j = 1,3
- z_prime(j) = -uz(j,i-1)
+ z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i)))
enddo
c write (2,*) "i",i
c write (2,*) "x_prime",(x_prime(j),j=1,3)
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 + dsign(1.0,dfloat(itype(i)))
- & *z_prime(j)*dc_norm(j,i+nres)
+ zz = zz + z_prime(j)*dc_norm(j,i+nres)
enddo
xxtab(i)=xx
Cc diagnostics - remove later
xx1 = dcos(alph(2))
yy1 = dsin(alph(2))*dcos(omeg(2))
- zz1 = -dsign(1.0,itype(i))*dsin(alph(2))*dsin(omeg(2))
+ zz1 = -dsign(1.0,dfloat(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
c & dscp1,dscp2,sumene
c sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
escloc = escloc + sumene
-c write (2,*) "i",i," escloc",sumene,escloc
+c write (2,*) "i",i," escloc",sumene,escloc,it,itype(i)
+c & ,zz,xx,yy
+c#define DEBUG
#ifdef DEBUG
C
C This section to check the numerical derivatives of the energy of ith side
C
C Compute the gradient of esc
C
+c zz=zz*dsign(1.0,dfloat(itype(i)))
pom_s1=(1.0d0+x(63))/(0.1d0 + dscp1)**2
pom_s16=6*(1.0d0+x(64))/(0.1d0 + dscp1**6)**2
pom_s2=(1.0d0+x(65))/(0.1d0 + dscp2)**2
& +(sumene2x+sumene4x*cost2tab(i+1))*(s2+s2_6)
& +(pom1+pom2)*pom_dx
#ifdef DEBUG
- write(2,*), "de_dxx = ", de_dxx,de_dxx_num
+ write(2,*), "de_dxx = ", de_dxx,de_dxx_num,itype(i)
#endif
C
sumene1y=x(3) + 2*x(6)*yy + x(9)*xx + x(10)*zz
& +(sumene2y+sumene4y*cost2tab(i+1))*(s2+s2_6)
& +(pom1-pom2)*pom_dy
#ifdef DEBUG
- write(2,*), "de_dyy = ", de_dyy,de_dyy_num
+ write(2,*), "de_dyy = ", de_dyy,de_dyy_num,itype(i)
#endif
C
de_dzz =(x(24) +2*x(27)*zz +x(28)*xx +x(30)*yy
& +x(60)*xx*yy)*cost2tab(i+1)*(s2+s2_6)
& + ( x(14) + 2*x(17)*zz+ x(18)*xx + x(20)*yy)*(s2+s2_6)
#ifdef DEBUG
- write(2,*), "de_dzz = ", de_dzz,de_dzz_num
+ write(2,*), "de_dzz = ", de_dzz,de_dzz_num,itype(i)
#endif
C
de_dt = 0.5d0*sumene3*cost2tab(i+1)*(s1+s1_6)
& -0.5d0*sumene4*sint2tab(i+1)*(s2+s2_6)
& +pom1*pom_dt1+pom2*pom_dt2
#ifdef DEBUG
- write(2,*), "de_dt = ", de_dt,de_dt_num
+ write(2,*), "de_dt = ", de_dt,de_dt_num,itype(i)
#endif
+c#undef DEBUG
c
C
cossc=scalar(dc_norm(1,i),dc_norm(1,i+nres))
dZZ_Ci1(k)=0.0d0
dZZ_Ci(k)=0.0d0
do j=1,3
- dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres)
- dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres)
+ dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)
+ & *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
+ dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)
+ & *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
enddo
dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
dYY_XYZ(k)=vbld_inv(i+nres)*(y_prime(k)-yy*dC_norm(k,i+nres))
- dZZ_XYZ(k)=vbld_inv(i+nres)*(z_prime(k)-zz*dC_norm(k,i+nres))
+ dZZ_XYZ(k)=vbld_inv(i+nres)*
+ & (z_prime(k)-zz*dC_norm(k,i+nres))
c
dt_dCi(k) = -dt_dCi(k)/sinttab(i+1)
dt_dCi1(k)= -dt_dCi1(k)/sinttab(i+1)
difi=phii-phi0(i)
if (difi.gt.drange(i)) then
difi=difi-drange(i)
- edihcnstr=edihcnstr+0.25d0*ftors*difi**4
- gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
else if (difi.lt.-drange(i)) then
difi=difi+drange(i)
- edihcnstr=edihcnstr+0.25d0*ftors*difi**4
- gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+ edihcnstr=edihcnstr+0.25d0*ftors(i)**difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
endif
! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
difi=pinorm(phii-phi0(i))
if (difi.gt.drange(i)) then
difi=difi-drange(i)
- edihcnstr=edihcnstr+0.25d0*ftors*difi**4
- gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
else if (difi.lt.-drange(i)) then
difi=difi+drange(i)
- edihcnstr=edihcnstr+0.25d0*ftors*difi**4
- gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
else
difi=0.0
endif
-cd write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii,
-cd & rad2deg*phi0(i), rad2deg*drange(i),
-cd & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
+ if (energy_dec) then
+ write (iout,'(a6,2i5,4f8.3,2e14.5)') "edihc",
+ & i,itori,rad2deg*phii,
+ & rad2deg*phi0(i), rad2deg*drange(i),
+ & rad2deg*difi,0.25d0*ftors(i)*difi**4,gloc(itori-3,icg)
+ endif
enddo
cd write (iout,*) 'edihcnstr',edihcnstr
return
cold ghalf=0.5d0*eel5*eij*gacont_hbr(ll,kk,k)
cgrad ghalf=0.5d0*ggg2(ll)
cd ghalf=0.0d0
- gradcorr5(ll,k)=gradcorr5(ll,k)+ghalf+ekont*derx(ll,2,2)
+ gradcorr5(ll,k)=gradcorr5(ll,k)+ekont*derx(ll,2,2)
gradcorr5(ll,k+1)=gradcorr5(ll,k+1)+ekont*derx(ll,3,2)
- gradcorr5(ll,l)=gradcorr5(ll,l)+ghalf+ekont*derx(ll,4,2)
+ gradcorr5(ll,l)=gradcorr5(ll,l)+ekont*derx(ll,4,2)
gradcorr5(ll,l1)=gradcorr5(ll,l1)+ekont*derx(ll,5,2)
gradcorr5_long(ll,l)=gradcorr5_long(ll,l)+gradcorr5kl
gradcorr5_long(ll,k)=gradcorr5_long(ll,k)-gradcorr5kl
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
C o o C
C \ /l\ /j\ / C
C \ / \ / \ / C
-C o| o | | o |o 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
+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,
double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2)
logical swap
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C C
+C C
C Parallel Antiparallel C
C C
-C o o C
+C o o C
C /l\ / \ /j\ C
C / \ / \ / \ C
C /| o |o o| o |\ C
& auxvec1(2),auxmat1(2,2)
logical swap
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C C
+C C
C Parallel Antiparallel C
C C
C o o C
C / \ / \ / \ C
C /| o |o o| o |\ C
C \ j|/k\| \ |/k\|l C
-C \ / \ \ / \ C
+C \ / \ \ / \ C
C o \ o \ C
C i i C
-C C
+C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C 4/7/01 AL Component s1 was removed, because it pertains to the respective
projects / unres.git / blobdiff