etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2+welec*fact(1)*ees
& +wvdwpp*evdw1
& +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+ & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
& +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
& +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2
& +welec*fact(1)*(ees+evdw1)
& +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+ & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
& +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
& +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
energia(19)=esccor
energia(20)=edihcnstr
energia(21)=evdw_t
+c if (dyn_ss) call dyn_set_nss
c detecting NaNQ
#ifdef ISNAN
#ifdef AIX
include 'COMMON.ENEPS'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
+ include 'COMMON.SBRIDGE'
logical lprn
common /srutu/icall
integer icant
C
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
+C in case of diagnostics write (iout,*) "TU SZUKAJ",i,j,dyn_ss_mask(i),dyn_ss_mask(j)
+C /06/28/2013 Adasko: In case of dyn_ss - dynamic disulfide bond
+C formation no electrostatic interactions should be calculated. If it
+C would be allowed NaN would appear
+ IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+C /06/28/2013 Adasko: dyn_ss_mask is logical statement wheather this Cys
+C residue can or cannot form disulfide bond. There is still bug allowing
+C Cys...Cys...Cys bond formation
+ call dyn_ssbond_ene(i,j,evdwij)
+C /06/28/2013 Adasko: dyn_ssbond_ene is dynamic SS bond foration energy
+C function in ssMD.F
+ evdw=evdw+evdwij
+c if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+c & 'evdw',i,j,evdwij,' ss'
+ ELSE
ind=ind+1
itypj=itype(j)
dscj_inv=vbld_inv(j+nres)
c write (iout,*) "i",i," j",j," itypi",itypi," itypj",itypj,
c & " ij",ij," eneps",aux*e1/dabs(eps(itypi,itypj)),
c & aux*e2/eps(itypi,itypj)
+c write (iout,'(a6,2i5,0pf7.3)') 'evdw',i,j,evdwij
if (lprn) then
sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
C Calculate angular part of the gradient.
call sc_grad
endif
+ ENDIF ! dyn_ss
enddo ! j
enddo ! iint
enddo ! i
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 (.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
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
+c & +akct*deltad*deltat12
& +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,
-c & " deltat12",deltat12," eij",eij
+ write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth,
+ & " akct",akct," deltad",deltad," deltat",deltat1,deltat2,
+ & " deltat12",deltat12," eij",eij,"ebr",ebr
ed=2*akcm*deltad+akct*deltat12
pom1=akct*deltad
pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
include 'COMMON.FFIELD'
include 'COMMON.CONTROL'
double precision u(3),ud(3)
+ logical :: lprn=.false.
estr=0.0d0
do i=nnt+1,nct
diff = vbld(i)-vbldp0
nbi=nbondterm(iti)
if (nbi.eq.1) then
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 (lprn)
+ & write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
+ & AKSC(1,iti),AKSC(1,iti)*diff*diff
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)
usum=usum+uprod1
usumsqder=usumsqder+ud(j)*uprod2
enddo
-c write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti),
-c & AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi)
+ if (lprn)
+ & write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti),
+ & AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi)
estr=estr+uprod/usum
do j=1,3
gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres)
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
escloc = escloc + sumene
c write (2,*) "escloc",escloc
if (.not. calc_grad) goto 1
-#ifdef DEBUG
+
+#ifdef DEBUG2
C
C This section to check the numerical derivatives of the energy of ith side
C chain in xx, yy, zz, and theta. Use the -DDEBUG compiler option or insert
itori=itortyp(itype(i-2))
itori1=itortyp(itype(i-1))
itori2=itortyp(itype(i))
-c iblock=1
-c if (iabs(itype(i+1)).eq.20) iblock=2
phii=phi(i)
phii1=phi(i+1)
gloci1=0.0D0
gloci2=0.0D0
C Regular cosine and sine terms
-c c do j=1,ntermd_1(itori,itori1,itori2,iblock)
-c v1cij=v1c(1,j,itori,itori1,itori2,iblock)
-c v1sij=v1s(1,j,itori,itori1,itori2,iblock)
-c v2cij=v1c(2,j,itori,itori1,itori2,iblock)
-c v2sij=v1s(2,j,itori,itori1,itori2,iblock)
- do j=1,ntermd_1(itori,itori1,itori2)
+ do j=1,ntermd_1(itori,itori1,itori2)
v1cij=v1c(1,j,itori,itori1,itori2)
v1sij=v1s(1,j,itori,itori1,itori2)
v2cij=v1c(2,j,itori,itori1,itori2)
v2sij=v1s(2,j,itori,itori1,itori2)
-
cosphi1=dcos(j*phii)
sinphi1=dsin(j*phii)
cosphi2=dcos(j*phii1)
gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
enddo
do k=2,ntermd_2(itori,itori1,itori2)
-c do k=2,ntermd_2(itori,itori1,itori2,iblock)
do l=1,k-1
-c v1cdij = v2c(k,l,itori,itori1,itori2,iblock)
-c v2cdij = v2c(l,k,itori,itori1,itori2,iblock)
-c v1sdij = v2s(k,l,itori,itori1,itori2,iblock)
-c v2sdij = v2s(l,k,itori,itori1,itori2,iblock)
v1cdij = v2c(k,l,itori,itori1,itori2)
v2cdij = v2c(l,k,itori,itori1,itori2)
v1sdij = v2s(k,l,itori,itori1,itori2)
C Set lprn=.true. for debugging
lprn=.false.
c lprn=.true.
-c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
+c write (iout,*) "EBACK_SC_COR",itau_start,itau_end,nterm_sccor
esccor=0.0D0
do i=itau_start,itau_end
esccor_ii=0.0D0
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