cmc
cmc Sep-06: egb takes care of dynamic ss bonds too
cmc
- if (dyn_ss) call dyn_set_nss
+c if (dyn_ss) call dyn_set_nss
c print *,"Processor",myrank," computed USCSC"
#ifdef TIMING
energia(23)=evdw_m
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
#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
& edihcnstr,ebr*nss,
& Uconst,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)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'UCONST= ',1pE16.6,' (Constraint energy)'/
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 Loop over i,i+2 and i,i+3 pairs of the peptide groups
C
do i=iturn3_start,iturn3_end
+C if (itype(i).eq.21 .or. itype(i+1).eq.21
+C & .or. itype(i+2).eq.21 .or. itype(i+3).eq.21.or.itype(i+4).eq.21)
+C & cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
num_cont_hb(i)=num_conti
enddo
do i=iturn4_start,iturn4_end
+C if (itype(i).eq.21 .or. itype(i+1).eq.21
+C & .or. itype(i+2).eq.21 .or. itype(i+3).eq.21.or.itype(i+4).eq.21
+C & .or. itype(i+5).eq.21)
+C & cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
c
do i=iatel_s,iatel_e
+C if (itype(i).eq.21 .or. itype(i+1).eq.21
+C &.or.itype(i+2)) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
+C if (itype(j).eq.21 .or. itype(j+1).eq.21
+C &.or.itype(j+2)) cycle
call eelecij(i,j,ees,evdw1,eel_loc)
enddo ! j
num_cont_hb(i)=num_conti
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,
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)
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) write (iout,*)
+ & "estr sc",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)
logical lprn /.false./, lprn1 /.false./
etheta=0.0D0
do i=ithet_start,ithet_end
+ if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+ &(itype(i).eq.ntyp1)) cycle
dethetai=0.0d0
dephii=0.0d0
dephii1=0.0d0
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
- if (i.gt.3) then
+C if (i.gt.3) then
+ 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
+ ityp1=ithetyp(itype(i-2))
do k=1,nsingle
cosph1(k)=0.0d0
sinph1(k)=0.0d0
enddo
endif
- if (i.lt.nres) 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
enddo
enddo
10 continue
- if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)')
- & i,theta(i)*rad2deg,phii*rad2deg,
+ if (lprn1) write (iout,'(a4,i2,3f8.1,9h ethetai ,f10.5)')
+ & 'ebe', i,theta(i)*rad2deg,phii*rad2deg,
& phii1*rad2deg,ethetai
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)=wang*dethetai
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))
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
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)
c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
esccor=0.0D0
do i=itau_start,itau_end
+C do i=42,42
esccor_ii=0.0D0
+ if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
isccori=isccortyp(itype(i-2))
isccori1=isccortyp(itype(i-1))
phii=phi(i)
+
cccc Added 9 May 2012
cc Tauangle is torsional engle depending on the value of first digit
c(see comment below)
cc Omicron is flat angle depending on the value of first digit
c(see comment below)
-
+C print *,i,tauangle(1,i)
- do intertyp=1,3 !intertyp
+c do intertyp=1,3 !intertyp
+ do intertyp=2,2 !intertyp
cc Added 09 May 2012 (Adasko)
cc Intertyp means interaction type of backbone mainchain correlation:
c 1 = SC...Ca...Ca...Ca
esccor=esccor+v1ij*cosphi+v2ij*sinphi
gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
enddo
+C print *,i,tauangle(1,i),gloci
gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
c write (iout,*) "WTF",intertyp,i,itype(i),v1ij*cosphi+v2ij*sinphi
c &gloc_sc(intertyp,i-3,icg)
enddo !intertyp
enddo
c do i=1,nres
-c write (iout,*) "W@T@F", gloc_sc(1,i,icg),gloc(i,icg)
+c write (iout,*) "W@T@F", gloc_sc(1,i,icg),gloc_sc(2,i,icg),
+c & gloc_sc(3,i,icg)
c enddo
return
end