& +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
& +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
- write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)')
- & i,(gvdwx(j,i),j=1,3),(gvdwc(j,i),j=1,3)
+ write (iout,'(i3,3f10.5,5x,3f10.5,5x,3f10.5,5x,3f10.5)')
+ & i,(gvdwx(j,i),j=1,3),(gvdwcT(j,i),j=1,3),(gvdwc(j,i),j=1,3),
+ & (gvdwcT(j,i),j=1,3)
+#ifdef TSCSC
+ do i=1,nct
+ do j=1,3
+ gradbufc(j,i)=wsc*gvdwc(j,i)+wsc*wscT*gvdwcT(j,i)+
+ & wscp*(gvdwc_scp(j,i)+gvdwc_scpp(j,i))+
+ & welec*gelc_long(j,i)+wvdwpp*gvdwpp(j,i)+
+ & wel_loc*gel_loc_long(j,i)+
+ & wcorr*gradcorr_long(j,i)+
+ & wcorr5*gradcorr5_long(j,i)+
+ & wcorr6*gradcorr6_long(j,i)+
+ & wturn6*gcorr6_turn_long(j,i)+
+ & wstrain*ghpbc(j,i)
+ enddo
+ enddo
+#else
+c call MPI_AllReduce(gradbufc(1,1),gradbufc_sum(1,1),3*nres,
+c & MPI_DOUBLE_PRECISION,MPI_SUM,FG_COMM,IERR)
+c time_reduce=time_reduce+MPI_Wtime()-time00
+#ifdef DEBUG
+c write (iout,*) "gradbufc_sum after allreduce"
+c do i=1,nres
+c write (iout,'(i3,3f10.5)') i,(gradbufc_sum(j,i),j=1,3)
+c enddo
+c call flush(iout)
+#endif
+c
+c Obsolete and inefficient code; we can make the effort O(n) and, therefore,
+c do not parallelize this part.
+c
+c do i=igrad_start,igrad_end
+c do j=jgrad_start(i),jgrad_end(i)
+c do k=1,3
+c gradbufc(k,i)=gradbufc(k,i)+gradbufc_sum(k,j)
+c enddo
+c enddo
+c enddo
gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
& wel_loc*gel_loc(j,i)+
& 0.5d0*(wscp*gvdwc_scpp(j,i)+
gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
& wel_loc*gel_loc(j,i)+
& 0.5d0*(wscp*gvdwc_scpp(j,i)+
gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
& wbond*gradbx(j,i)+
& wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
& wsccor*gsccorx(j,i)
& +wscloc*gsclocx(j,i)
gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
& wbond*gradbx(j,i)+
& wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
& wsccor*gsccorx(j,i)
& +wscloc*gsclocx(j,i)
& MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres,
& MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
& MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres,
& MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
gvdwc_scp_norm=dsqrt(scalar(gvdwc_scp(1,i),gvdwc_scp(1,i)))
if (gvdwc_scp_norm.gt.gvdwc_scp_max)
& gvdwc_scp_max=gvdwc_scp_norm
gvdwc_scp_norm=dsqrt(scalar(gvdwc_scp(1,i),gvdwc_scp(1,i)))
if (gvdwc_scp_norm.gt.gvdwc_scp_max)
& gvdwc_scp_max=gvdwc_scp_norm
if (gscloc_norm.gt.gscloc_max) gscloc_max=gscloc_norm
gvdwx_norm=dsqrt(scalar(gvdwx(1,i),gvdwx(1,i)))
if (gvdwx_norm.gt.gvdwx_max) gvdwx_max=gvdwx_norm
if (gscloc_norm.gt.gscloc_max) gscloc_max=gscloc_norm
gvdwx_norm=dsqrt(scalar(gvdwx(1,i),gvdwx(1,i)))
if (gvdwx_norm.gt.gvdwx_max) gvdwx_max=gvdwx_norm
gradx_scp_norm=dsqrt(scalar(gradx_scp(1,i),gradx_scp(1,i)))
if (gradx_scp_norm.gt.gradx_scp_max)
& gradx_scp_max=gradx_scp_norm
gradx_scp_norm=dsqrt(scalar(gradx_scp(1,i),gradx_scp(1,i)))
if (gradx_scp_norm.gt.gradx_scp_max)
& gradx_scp_max=gradx_scp_norm
- & '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)'/
- & 'EHBP= ',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,
- & '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)'/
& '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
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
-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)
+ else if (ii.gt.nres .and. jj.gt.nres) then
+c Restraints from contact prediction
+ dd=dist(ii,jj)
+ 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
+ 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
- dd=dist(ii,jj)
- rdis=dd-dhpb(i)
+ dd=dist(ii,jj)
+ 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)
cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
cd & ' waga=',waga,' fac=',fac
cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
cd & ' waga=',waga,' fac=',fac
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).
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).
& +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,
& +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,
- 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,
& restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
& (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6)
gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
& restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
& (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6)
gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
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)
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)
do j=1,ntermd_1(itori,itori1,itori2)
v1cij=v1c(1,j,itori,itori1,itori2)
v1sij=v1s(1,j,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)
+
+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)
+
+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
+c 2 = Ca...Ca...Ca...SC
+c 3 = SC...Ca...Ca...SCi
- do j=1,nterm_sccor
- v1ij=v1sccor(j,itori,itori1)
- v2ij=v2sccor(j,itori,itori1)
- cosphi=dcos(j*phii)
- sinphi=dsin(j*phii)
+ if (((intertyp.eq.3).and.((itype(i-2).eq.10).or.
+ & (itype(i-1).eq.10).or.(itype(i-2).eq.21).or.
+ & (itype(i-1).eq.21)))
+ & .or. ((intertyp.eq.1).and.((itype(i-2).eq.10)
+ & .or.(itype(i-2).eq.21)))
+ & .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or.
+ & (itype(i-1).eq.21)))) cycle
+ if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.21)) cycle
+ if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.21))
+ & cycle
+ do j=1,nterm_sccor(isccori,isccori1)
+ v1ij=v1sccor(j,intertyp,isccori,isccori1)
+ v2ij=v2sccor(j,intertyp,isccori,isccori1)
+ cosphi=dcos(j*tauangle(intertyp,i))
+ sinphi=dsin(j*tauangle(intertyp,i))
+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)
if (lprn)
& write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
& restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
if (lprn)
& write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
& restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
- & (v1sccor(j,itori,itori1),j=1,6),(v2sccor(j,itori,itori1),j=1,6)
+ & (v1sccor(j,intertyp,itori,itori1),j=1,6)
+ & ,(v2sccor(j,intertyp,itori,itori1),j=1,6)
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,
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
double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2)
logical swap
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC