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
#endif
#ifdef MPI
include 'mpif.h'
+#endif
double precision gradbufc(3,maxres),gradbufx(3,maxres),
& glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres)
-#endif
include 'COMMON.SETUP'
include 'COMMON.IOUNITS'
include 'COMMON.FFIELD'
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
C Compute the virtual-bond-torsional-angle dependent quantities needed
C to calculate the el-loc multibody terms of various order.
C
- write(iout,*) 'nphi=',nphi,nres
+c write(iout,*) 'nphi=',nphi,nres
#ifdef PARMAT
do i=ivec_start+2,ivec_end+2
#else
else
iti1=ntortyp+1
endif
- write(iout,*),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 write(iout,*),i
+ b1(1,i-2)=bnew1(1,1,iti)*dsin(theta(i-1)/2.0)
+ & +bnew1(2,1,iti)*dsin(theta(i-1))
+ & +bnew1(3,1,iti)*dcos(theta(i-1)/2.0)
+ 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)*sin(alpha(i))*cos(beta(i))
c &*(cos(theta(i)/2.0)
- 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)
+ b2(1,i-2)=bnew2(1,1,iti)*dsin(theta(i-1)/2.0)
+ & +bnew2(2,1,iti)*dsin(theta(i-1))
+ & +bnew2(3,1,iti)*dcos(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)*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
+ 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
c if (ggb1(1,i).eq.0.0d0) then
c write(iout,*) 'i=',i,ggb1(1,i),
c &bnew1(1,1,iti)*cos(theta(i)/2.0)/2.0,
gtb1(2,i-2)=0.0
b2(2,i-2)=bnew2(1,2,iti)
gtb2(2,i-2)=0.0
-c EE(1,1,iti)=0.0d0
+ EE(1,1,i-2)=eenew(1,iti)*dcos(theta(i-1))
+ EE(1,2,i-2)=eeold(1,2,iti)
+ EE(2,1,i-2)=eeold(2,1,iti)
+ EE(2,2,i-2)=eeold(2,2,iti)
+ gtEE(1,1,i-2)=-eenew(1,iti)*dsin(theta(i-1))
+ gtEE(1,2,i-2)=0.0d0
+ gtEE(2,2,i-2)=0.0d0
+ gtEE(2,1,i-2)=0.0d0
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)
b1tilde(2,i-2)=-b1(2,i-2)
b2tilde(1,i-2)=b2(1,i-2)
b2tilde(2,i-2)=-b2(2,i-2)
- write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2)
- write (iout,*) 'theta=', theta(i-1)
+c write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2)
+c write(iout,*) 'b1=',b1(1,i-2)
+c write (iout,*) 'theta=', theta(i-1)
enddo
#ifdef PARMAT
do i=ivec_start+2,ivec_end+2
Ug2der(2,2,i-2)=0.0d0
endif
c if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then
-#ifndef NEWCORR
if (i.gt. nnt+2 .and. i.lt.nct+2) then
iti = itortyp(itype(i-2))
else
else
iti1=ntortyp+1
endif
-#endif
cd write (iout,*) '*******i',i,' iti1',iti
cd write (iout,*) 'b1',b1(:,iti)
cd write (iout,*) 'b2',b2(:,iti)
call matvec2(Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2))
c write (iout,*) Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2),"chuj"
#endif
- call matmat2(EE(1,1,iti),Ug(1,1,i-2),EUg(1,1,i-2))
+c write(iout,*) "co jest kurwa", iti, EE(1,1,iti),EE(2,1,iti),
+c & EE(1,2,iti),EE(2,2,iti)
+ call matmat2(EE(1,1,i-2),Ug(1,1,i-2),EUg(1,1,i-2))
+ call matmat2(gtEE(1,1,i-2),Ug(1,1,i-2),gtEUg(1,1,i-2))
+c write(iout,*) "Macierz EUG",
+c & eug(1,1,i-2),eug(1,2,i-2),eug(2,1,i-2),
+c & eug(2,2,i-2)
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0)
& then
call matmat2(CC(1,1,iti),Ug(1,1,i-2),CUg(1,1,i-2))
enddo
endif
call matvec2(Ugder(1,1,i-2),b2(1,i-2),Ub2der(1,i-2))
- call matmat2(EE(1,1,iti),Ugder(1,1,i-2),EUgder(1,1,i-2))
+ call matmat2(EE(1,1,i-2),Ugder(1,1,i-2),EUgder(1,1,i-2))
do k=1,2
muder(k,i-2)=Ub2der(k,i-2)
enddo
do k=1,2
mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1)
enddo
-cd write (iout,*) 'mu ',mu(:,i-2)
+c write (iout,*) 'mu ',mu(:,i-2),i-2
cd write (iout,*) 'mu1',mu1(:,i-2)
cd write (iout,*) 'mu2',mu2(:,i-2)
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0)
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
num_conti=0
-c TU ZLE
-c call eelecij(i,i+2,ees,evdw1,eel_loc)
+ call eelecij(i,i+2,ees,evdw1,eel_loc)
if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
num_cont_hb(i)=num_conti
enddo
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
num_conti=num_cont_hb(i)
-c TU ZLE
-c call eelecij(i,i+3,ees,evdw1,eel_loc)
+c write(iout,*) "JESTEM W PETLI"
+ call eelecij(i,i+3,ees,evdw1,eel_loc)
if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
& call eturn4(i,eello_turn4)
num_cont_hb(i)=num_conti
c
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
c
-c do i=iatel_s,iatel_e
- do i=4,5
+ do i=iatel_s,iatel_e
+c do i=7,7
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
zmedi=c(3,i)+0.5d0*dzi
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
-c do j=ielstart(i),ielend(i)
- do j=8,9
- write (iout,*) 'tu wchodze',i,j,itype(i),itype(j)
+ do j=ielstart(i),ielend(i)
+c do j=13,13
+c write (iout,*) 'tu wchodze',i,j,itype(i),itype(j)
if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle
call eelecij(i,j,ees,evdw1,eel_loc)
enddo ! j
do l=1,2
kkk=kkk+1
muij(kkk)=mu(k,i)*mu(l,j)
+c write(iout,*) 'mumu=', mu(k,i),mu(l,j),i,j,k,l
#ifdef NEWCORR
- gmuij1(kkk)=gtb1(k,i)*mu(l,j)
- write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(k,i),k,i
- gmuij2(kkk)=gUb2(k,i-1)*mu(l,j)
- gmuji1(kkk)=mu(k,i)*gtb1(l,j)
- write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(l,j),l,j
- gmuji2(kkk)=mu(k,i)*gUb2(l,j-1)
+ gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
+c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
+ gmuij2(kkk)=gUb2(k,i)*mu(l,j)
+ gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
+c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
+ gmuji2(kkk)=mu(k,i)*gUb2(l,j)
#endif
enddo
enddo
C Contribution to the local-electrostatic energy coming from the i-j pair
eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
& +a33*muij(4)
+c write(iout,*) 'muije=',muij(1),muij(2),muij(3),muij(4)
C Calculate patrial derivative for theta angle
#ifdef NEWCORR
geel_loc_ij=a22*gmuij1(1)
& +a23*gmuij1(2)
& +a32*gmuij1(3)
& +a33*gmuij1(4)
- write(iout,*) "derivative over thatai"
- write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
- & a33*gmuij1(4)
+c write(iout,*) "derivative over thatai"
+c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
+c & a33*gmuij1(4)
gloc(nphi+i,icg)=gloc(nphi+i,icg)+
& geel_loc_ij*wel_loc
- write(iout,*) "derivative over thatai-1"
- write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
- & a33*gmuij2(4)
- geel_loc_ij=a22*gmuij2(1)+a23*gmuij2(2)+a32*gmuij2(3)
+c write(iout,*) "derivative over thatai-1"
+c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
+c & a33*gmuij2(4)
+ geel_loc_ij=
+ & a22*gmuij2(1)
+ & +a23*gmuij2(2)
+ & +a32*gmuij2(3)
& +a33*gmuij2(4)
gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
& geel_loc_ij*wel_loc
- geel_loc_ji=a22*gmuji1(1)+a23*gmuji1(2)+a32*gmuji1(3)
+c Derivative over j residue
+ geel_loc_ji=a22*gmuji1(1)
+ & +a23*gmuji1(2)
+ & +a32*gmuji1(3)
& +a33*gmuji1(4)
- write(iout,*) "derivative over thataj"
- write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
- & a33*gmuji1(4)
+c write(iout,*) "derivative over thataj"
+c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
+c & a33*gmuji1(4)
- gloc(nphi+j,icg)=gloc(nphi+j,icg)+
+ gloc(nphi+j,icg)=gloc(nphi+j,icg)+
& geel_loc_ji*wel_loc
- geel_loc_ji=a22*gmuji2(1)+a23*gmuji2(2)+a32*gmuji2(3)
+ geel_loc_ji=
+ & +a22*gmuji2(1)
+ & +a23*gmuji2(2)
+ & +a32*gmuji2(3)
& +a33*gmuji2(4)
- write(iout,*) "derivative over thataj-1"
- write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
- & a33*gmuji2(4)
+c write(iout,*) "derivative over thataj-1"
+c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
+c & a33*gmuji2(4)
gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
& geel_loc_ji*wel_loc
#endif
dimension ggg(3)
double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2),
& e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2),
- & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2)
+ & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),gpizda1(2,2),
+ & gpizda2(2,2),auxgmat1(2,2),auxgmatt1(2,2),
+ & auxgmat2(2,2),auxgmatt2(2,2)
double precision agg(3,4),aggi(3,4),aggi1(3,4),
& aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2)
common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
cd call checkint_turn3(i,a_temp,eello_turn3_num)
call matmat2(EUg(1,1,i+1),EUg(1,1,i+2),auxmat(1,1))
+c auxalary matices for theta gradient
+c auxalary matrix for i+1 and constant i+2
+ call matmat2(gtEUg(1,1,i+1),EUg(1,1,i+2),auxgmat1(1,1))
+c auxalary matrix for i+2 and constant i+1
+ call matmat2(EUg(1,1,i+1),gtEUg(1,1,i+2),auxgmat2(1,1))
call transpose2(auxmat(1,1),auxmat1(1,1))
+ call transpose2(auxgmat1(1,1),auxgmatt1(1,1))
+ call transpose2(auxgmat2(1,1),auxgmatt2(1,1))
call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
+ call matmat2(a_temp(1,1),auxgmatt1(1,1),gpizda1(1,1))
+ call matmat2(a_temp(1,1),auxgmatt2(1,1),gpizda2(1,1))
eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2))
+C Derivatives in theta
+ gloc(nphi+i,icg)=gloc(nphi+i,icg)
+ & +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3
+ gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)
+ & +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3
+
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
& 'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2))
cd write (2,*) 'i,',i,' j',j,'eello_turn3',
dimension ggg(3)
double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2),
& e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2),
- & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),auxgvec(2)
+ & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),auxgvec(2),
+ & auxgEvec1(2),auxgEvec2(2),auxgEvec3(2),
+ & gte1t(2,2),gte2t(2,2),gte3t(2,2),
+ & gte1a(2,2),gtae3(2,2),gtae3e2(2,2), ae3gte2(2,2),
+ & gtEpizda1(2,2),gtEpizda2(2,2),gtEpizda3(2,2)
double precision agg(3,4),aggi(3,4),aggi1(3,4),
& aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2)
common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
cd call checkint_turn4(i,a_temp,eello_turn4_num)
c write (iout,*) "eturn4 i",i," j",j," j1",j1," j2",j2
+c write(iout,*)"WCHODZE W PROGRAM"
a_temp(1,1)=a22
a_temp(1,2)=a23
a_temp(2,1)=a32
call transpose2(EUg(1,1,i+1),e1t(1,1))
call transpose2(Eug(1,1,i+2),e2t(1,1))
call transpose2(Eug(1,1,i+3),e3t(1,1))
+C Ematrix derivative in theta
+ call transpose2(gtEUg(1,1,i+1),gte1t(1,1))
+ call transpose2(gtEug(1,1,i+2),gte2t(1,1))
+ call transpose2(gtEug(1,1,i+3),gte3t(1,1))
call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
+c eta1 in derivative theta
+ call matmat2(gte1t(1,1),a_temp(1,1),gte1a(1,1))
call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
+c auxgvec is derivative of Ub2 so i+3 theta
call matvec2(e1a(1,1),gUb2(1,i+3),auxgvec(1))
+c auxalary matrix of E i+1
+ call matvec2(gte1a(1,1),Ub2(1,i+3),auxgEvec1(1))
c s1=0.0
c gs1=0.0
s1=scalar2(b1(1,i+2),auxvec(1))
-c gs1=scalar2(gtb1(1,i+2),auxgvec(1))
+c derivative of theta i+2 with constant i+3
+ gs23=scalar2(gtb1(1,i+2),auxvec(1))
+c derivative of theta i+2 with constant i+2
+ gs32=scalar2(b1(1,i+2),auxgvec(1))
+c derivative of E matix in theta of i+1
+ gsE13=scalar2(b1(1,i+2),auxgEvec1(1))
+
call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
+c ea31 in derivative theta
+ call matmat2(a_temp(1,1),gte3t(1,1),gtae3(1,1))
call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
+c auxilary matrix auxgvec of Ub2 with constant E matirx
call matvec2(ae3(1,1),gUb2(1,i+2),auxgvec(1))
+c auxilary matrix auxgEvec1 of E matix with Ub2 constant
+ call matvec2(gtae3(1,1),Ub2(1,i+2),auxgEvec3(1))
+
c s2=0.0
c gs2=0.0
s2=scalar2(b1(1,i+1),auxvec(1))
-c gs2=scalar2(gtb1(1,i+1),auxgvec(1))
-c write(iout,*) gs1,gs2,'i=',i,auxgvec(1),gUb2(1,i+2),ggb1(1,i+2),
-c & ggb1(1,i+1)
+c derivative of theta i+1 with constant i+3
+ gs13=scalar2(gtb1(1,i+1),auxvec(1))
+c derivative of theta i+2 with constant i+1
+ gs21=scalar2(b1(1,i+1),auxgvec(1))
+c derivative of theta i+3 with constant i+1
+ gsE31=scalar2(b1(1,i+1),auxgEvec3(1))
+c write(iout,*) gs1,gs2,'i=',i,auxgvec(1),gUb2(1,i+2),gtb1(1,i+2),
+c & gtb1(1,i+1)
call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
+c two derivatives over diffetent matrices
+c gtae3e2 is derivative over i+3
+ call matmat2(gtae3(1,1),e2t(1,1),gtae3e2(1,1))
+c ae3gte2 is derivative over i+2
+ call matmat2(ae3(1,1),gte2t(1,1),ae3gte2(1,1))
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
+c three possible derivative over theta E matices
+c i+1
+ call matmat2(ae3e2(1,1),gte1t(1,1),gtEpizda1(1,1))
+c i+2
+ call matmat2(ae3gte2(1,1),e1t(1,1),gtEpizda2(1,1))
+c i+3
+ call matmat2(gtae3e2(1,1),e1t(1,1),gtEpizda3(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
+
+ gsEE1=0.5d0*(gtEpizda1(1,1)+gtEpizda1(2,2))
+ gsEE2=0.5d0*(gtEpizda2(1,1)+gtEpizda2(2,2))
+ gsEE3=0.5d0*(gtEpizda3(1,1)+gtEpizda3(2,2))
+
eello_turn4=eello_turn4-(s1+s2+s3)
#ifdef NEWCORR
-c geel_loc_ij=-(gs1+gs2)
-c gloc(nphi+i,icg)=gloc(nphi+i,icg)-
-c & gs1
+ gloc(nphi+i,icg)=gloc(nphi+i,icg)
+ & -(gs13+gsE13+gsEE1)*wturn4
+ gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg)
+ & -(gs23+gs21+gsEE2)*wturn4
+ gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg)
+ & -(gs32+gsE31+gsEE3)*wturn4
c gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)-
c & gs2
#endif
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
& 'eturn4',i,j,-(s1+s2+s3)
-cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
-cd & ' eello_turn4_num',8*eello_turn4_num
+c write (iout,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
+c & ' eello_turn4_num',8*eello_turn4_num
C Derivatives in gamma(i)
call transpose2(EUgder(1,1,i+1),e1tder(1,1))
call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1))
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.
& iabs(itype(jjj)).eq.1) then
+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
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
C target distance.
- dd=dist(ii,jj)
- rdis=dd-dhpb(i)
+ dd=dist(ii,jj)
+ 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
C Evaluate gradient.
C
- fac=waga*rdis/dd
+ 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
+ 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
& 'ebend',i,ethetai
if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
- gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)
+ gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg)
enddo
C Ufff.... We've done all this!!!
return
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)=wang*dethetai+gloc(nphi+i-2,icg)
enddo
return
end
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