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()
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
do i=1,4*nres
glocbuf(i)=gloc(i,icg)
enddo
+#define DEBUG
#ifdef DEBUG
write (iout,*) "gloc_sc before reduce"
do i=1,nres
enddo
enddo
#endif
+#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
#ifdef DEBUG
write (iout,*) "gloc_sc after reduce"
do i=1,nres
enddo
enddo
#endif
+#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
evdw=0.0D0
ccccc energy_dec=.false.
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
+ 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=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
iti1=ntortyp+1
endif
c write(iout,*),i
- b1(1,i-2)=bnew1(1,1,iti)*dsin(theta(i-1)/2.0d0)
- & +bnew1(2,1,iti)*dsin(theta(i-1))
- & +bnew1(3,1,iti)*dcos(theta(i-1)/2.0d0)
- gtb1(1,i-2)=bnew1(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0
- & +bnew1(2,1,iti)*dcos(theta(i-1))
- & -bnew1(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0
-c & +bnew1(3,1,iti)*dsin(alpha(i))*cos(beta(i))
+ b1(1,i-2)=bnew1(1,1,iti)*sin(theta(i-1)/2.0)
+ & +bnew1(2,1,iti)*sin(theta(i-1))
+ & +bnew1(3,1,iti)*cos(theta(i-1)/2.0)
+ gtb1(1,i-2)=bnew1(1,1,iti)*cos(theta(i-1)/2.0)/2.0
+ & +bnew1(2,1,iti)*cos(theta(i-1))
+ & -bnew1(3,1,iti)*sin(theta(i-1)/2.0)/2.0
+c & +bnew1(3,1,iti)*sin(alpha(i))*cos(beta(i))
c &*(cos(theta(i)/2.0)
- b2(1,i-2)=bnew2(1,1,iti)*dsin(theta(i-1)/2.0d0)
- & +bnew2(2,1,iti)*dsin(theta(i-1))
- & +bnew2(3,1,iti)*dcos(theta(i-1)/2.0d0)
-c & +bnew2(3,1,iti)*dsin(alpha(i))*dcos(beta(i))
+ b2(1,i-2)=bnew2(1,1,iti)*sin(theta(i-1)/2.0)
+ & +bnew2(2,1,iti)*sin(theta(i-1))
+ & +bnew2(3,1,iti)*cos(theta(i-1)/2.0)
+c & +bnew2(3,1,iti)*sin(alpha(i))*cos(beta(i))
c &*(cos(theta(i)/2.0)
- gtb2(1,i-2)=bnew2(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0
- & +bnew2(2,1,iti)*dcos(theta(i-1))
- & -bnew2(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0
+ gtb2(1,i-2)=bnew2(1,1,iti)*cos(theta(i-1)/2.0)/2.0
+ & +bnew2(2,1,iti)*cos(theta(i-1))
+ & -bnew2(3,1,iti)*sin(theta(i-1)/2.0)/2.0
c if (ggb1(1,i).eq.0.0d0) then
c write(iout,*) 'i=',i,ggb1(1,i),
-c &bnew1(1,1,iti)*dcos(theta(i)/2.0d0)/2.0d0,
-c &bnew1(2,1,iti)*dcos(theta(i)),
-c &bnew1(3,1,iti)*dsin(theta(i)/2.0d0)/2.0d0
+c &bnew1(1,1,iti)*cos(theta(i)/2.0)/2.0,
+c &bnew1(2,1,iti)*cos(theta(i)),
+c &bnew1(3,1,iti)*sin(theta(i)/2.0)/2.0
c endif
b1(2,i-2)=bnew1(1,2,iti)
gtb1(2,i-2)=0.0
c EE(2,2,iti)=0.0d0
c EE(1,2,iti)=0.5d0*eenew(1,iti)
c EE(2,1,iti)=0.5d0*eenew(1,iti)
-c b1(2,iti)=bnew1(1,2,iti)*dsin(alpha(i))*dsin(beta(i))
-c b2(2,iti)=bnew2(1,2,iti)*dsin(alpha(i))*dsin(beta(i))
+c b1(2,iti)=bnew1(1,2,iti)*sin(alpha(i))*sin(beta(i))
+c b2(2,iti)=bnew2(1,2,iti)*sin(alpha(i))*sin(beta(i))
b1tilde(1,i-2)=b1(1,i-2)
b1tilde(2,i-2)=-b1(2,i-2)
b2tilde(1,i-2)=b2(1,i-2)
do i=3,nres+1
#endif
#endif
- if (i.gt. nnt+2 .and. i.lt.nct+2) then
- iti = itortyp(itype(i-2))
- else
- iti=ntortyp+1
- endif
-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))
- else
- iti1=ntortyp+1
- endif
if (i .lt. nres+1) then
sin1=dsin(phi(i))
cos1=dcos(phi(i))
mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1)
enddo
#ifdef MUOUT
- write (iout,'(2hmu,i3,3f8.1,12f10.5)') i-2,rad2deg*theta(i-1),
+ write (iout,'(2hmu,i3,3f8.1,7f10.5)') i-2,rad2deg*theta(i-1),
& rad2deg*theta(i),rad2deg*phi(i),mu(1,i-2),mu(2,i-2),
& -b2(1,i-2),b2(2,i-2),b1(1,i-2),b1(2,i-2),
& dsqrt(b2(1,i-1)**2+b2(2,i-1)**2)
- & +dsqrt(b1(1,i-1)**2+b1(2,i-1)**2),
- & ((ee(l,k,i-2),l=1,2),k=1,2),eenew(1,itortyp(iti))
+ & +dsqrt(b1(1,i-1)**2+b1(2,i-1)**2)
#endif
cd write (iout,*) 'mu ',mu(:,i-2)
cd write (iout,*) 'mu1',mu1(:,i-2)
cd write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj
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.
- & iabs(itype(jjj)).eq.1) then
- call ssbond_ene(iii,jjj,eij)
+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
+c if (.not.dyn_ss .and. i.le.nss) then
+C 15/02/13 CC dynamic SSbond
+C if (.not.dyn_ss.and.
+C & ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+
+ 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
C Calculate the distance between the two points and its difference from the
logical lprn
C Set lprn=.true. for debugging
lprn=.false.
-c lprn=.true.
+c lprn=.true.
etors=0.0D0
do i=iphi_start,iphi_end
if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
C Subtract the constant term
etors=etors-v0(itori,itori1,iblock)
if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
- & 'etor',i,etors_ii
+ & 'etor',i,etors_ii-v0(itori,itori1,iblock)
if (lprn)
- & write (iout,'(2(a3,2x,i3,2x),2i3,f10.2,6f8.3/36x,6f8.3/)')
+ & 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,
- & rad2deg*phii,
- & (v1(j,itori,itori1,1),j=1,6),(v2(j,itori,itori1,1),j=1,6)
+ & (v1(j,itori,itori1,iblock),j=1,6),
+ & (v2(j,itori,itori1,iblock),j=1,6)
gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
enddo
projects / unres.git / blobdiff