energia(17)=estr
energia(20)=Uconst+Uconst_back
energia(21)=esccor
+c Here are the energies showed per procesor if the are more processors
+c per molecule then we sum it up in sum_energy subroutine
c print *," Processor",myrank," calls SUM_ENERGY"
call sum_energy(energia,.true.)
c print *," Processor",myrank," left SUM_ENERGY"
#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
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
#ifdef PARMAT
do i=ivec_start+2,ivec_end+2
#else
do i=3,nres+1
#endif
+#ifdef NEWCORR
+ 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
+ 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 & +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)
+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
+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,
+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
+ b2(2,i-2)=bnew2(1,2,iti)
+ gtb2(2,i-2)=0.0
+ 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)
+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)
+ 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)
+ enddo
+#ifdef PARMAT
+ do i=ivec_start+2,ivec_end+2
+#else
+ do i=3,nres+1
+#endif
+#endif
if (i .lt. nres+1) then
sin1=dsin(phi(i))
cos1=dcos(phi(i))
cd write (iout,*) 'Ug',Ug(:,:,i-2)
c if (i .gt. iatel_s+2) then
if (i .gt. nnt+2) then
- call matvec2(Ug(1,1,i-2),b2(1,iti),Ub2(1,i-2))
- call matmat2(EE(1,1,iti),Ug(1,1,i-2),EUg(1,1,i-2))
+ call matvec2(Ug(1,1,i-2),b2(1,i-2),Ub2(1,i-2))
+#ifdef NEWCORR
+ 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
+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
enddo
endif
- call matvec2(Ugder(1,1,i-2),b2(1,iti),Ub2der(1,i-2))
- call matmat2(EE(1,1,iti),Ugder(1,1,i-2),EUgder(1,1,i-2))
+ call matvec2(Ugder(1,1,i-2),b2(1,i-2),Ub2der(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
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))
+ if (itype(i-1).le.ntyp) then
+ iti1 = itortyp(itype(i-1))
+ else
+ iti1=ntortyp+1
+ endif
else
iti1=ntortyp+1
endif
do k=1,2
- mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1)
+ mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1)
enddo
cd write (iout,*) 'mu ',mu(:,i-2)
cd write (iout,*) 'mu1',mu1(:,i-2)
call matvec2(Ctilde(1,1,iti1),obrot_der(1,i-2),Ctobrder(1,i-2))
call matvec2(Dtilde(1,1,iti),obrot2_der(1,i-2),Dtobr2der(1,i-2))
C Vectors and matrices dependent on a single virtual-bond dihedral.
- call matvec2(DD(1,1,iti),b1tilde(1,iti1),auxvec(1))
+ call matvec2(DD(1,1,iti),b1tilde(1,i-1),auxvec(1))
call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2))
call matvec2(Ug2der(1,1,i-2),auxvec(1),Ug2Db1tder(1,i-2))
call matvec2(CC(1,1,iti1),Ub2(1,i-2),CUgb2(1,i-2))
dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
& erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
- & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
+ & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij(4)
common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
& dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
& num_conti,j1,j2
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
num_conti=num_cont_hb(i)
+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)
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 do i=7,7
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,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 write (iout,*) i,j,itype(i),itype(j)
+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
dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
& erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
- & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
+ & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
+ & gmuij2(4),gmuji2(4)
common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
& dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
& num_conti,j1,j2
cd & xmedi,ymedi,zmedi,xj,yj,zj
if (energy_dec) then
- write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
+ write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)')
+ &'evdw1',i,j,evdwij
+ &,iteli,itelj,aaa,evdw1
write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
endif
C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
C are computed for EVERY pair of non-contiguous peptide groups.
C
+
if (j.lt.nres-1) then
j1=j+1
j2=j-1
j2=j-2
endif
kkk=0
+ lll=0
do k=1,2
do l=1,2
kkk=kkk+1
muij(kkk)=mu(k,i)*mu(l,j)
+#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)
+#endif
enddo
enddo
cd write (iout,*) 'EELEC: i',i,' j',j
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 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)
+ 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)
+ & +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)
+ & +a33*gmuji1(4)
+ write(iout,*) "derivative over thataj"
+ write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
+ & a33*gmuji1(4)
+
+ 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)
+ & +a33*gmuji2(4)
+ write(iout,*) "derivative over thataj-1"
+ write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
+ & a33*gmuji2(4)
+ gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
+ & geel_loc_ji*wel_loc
+#endif
cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
& 'eelloc',i,j,eel_loc_ij
+c write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
eel_loc=eel_loc+eel_loc_ij
C Partial derivatives in virtual-bond dihedral angles gamma
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)
+ & 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))
- s1=scalar2(b1(1,iti2),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 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))
- s2=scalar2(b1(1,iti1),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 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
+ 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))
call matvec2(auxmat(1,1),Ub2(1,i+3),auxvec(1))
- s1=scalar2(b1(1,iti2),auxvec(1))
+ s1=scalar2(b1(1,i+2),auxvec(1))
call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3)
C Derivatives in gamma(i+1)
call transpose2(EUgder(1,1,i+2),e2tder(1,1))
call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1))
- s2=scalar2(b1(1,iti1),auxvec(1))
+ s2=scalar2(b1(1,i+1),auxvec(1))
call matmat2(ae3(1,1),e2tder(1,1),auxmat(1,1))
call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
C Derivatives in gamma(i+2)
call transpose2(EUgder(1,1,i+3),e3tder(1,1))
call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1))
- s1=scalar2(b1(1,iti2),auxvec(1))
+ s1=scalar2(b1(1,i+2),auxvec(1))
call matmat2(a_temp(1,1),e3tder(1,1),auxmat(1,1))
call matvec2(auxmat(1,1),Ub2(1,i+2),auxvec(1))
- s2=scalar2(b1(1,iti1),auxvec(1))
+ s2=scalar2(b1(1,i+1),auxvec(1))
call matmat2(auxmat(1,1),e2t(1,1),auxmat3(1,1))
call matmat2(auxmat3(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
a_temp(2,2)=agg(l,4)
call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
- s1=scalar2(b1(1,iti2),auxvec(1))
+ s1=scalar2(b1(1,i+2),auxvec(1))
call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
- s2=scalar2(b1(1,iti1),auxvec(1))
+ s2=scalar2(b1(1,i+1),auxvec(1))
call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
a_temp(2,2)=aggi(l,4)
call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
- s1=scalar2(b1(1,iti2),auxvec(1))
+ s1=scalar2(b1(1,i+2),auxvec(1))
call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
- s2=scalar2(b1(1,iti1),auxvec(1))
+ s2=scalar2(b1(1,i+1),auxvec(1))
call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
a_temp(2,2)=aggi1(l,4)
call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
- s1=scalar2(b1(1,iti2),auxvec(1))
+ s1=scalar2(b1(1,i+2),auxvec(1))
call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
- s2=scalar2(b1(1,iti1),auxvec(1))
+ s2=scalar2(b1(1,i+1),auxvec(1))
call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
a_temp(2,2)=aggj(l,4)
call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
- s1=scalar2(b1(1,iti2),auxvec(1))
+ s1=scalar2(b1(1,i+2),auxvec(1))
call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
- s2=scalar2(b1(1,iti1),auxvec(1))
+ s2=scalar2(b1(1,i+1),auxvec(1))
call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
a_temp(2,2)=aggj1(l,4)
call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
- s1=scalar2(b1(1,iti2),auxvec(1))
+ s1=scalar2(b1(1,i+2),auxvec(1))
call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
- s2=scalar2(b1(1,iti1),auxvec(1))
+ s2=scalar2(b1(1,i+1),auxvec(1))
call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
cosphi=om12-om1*om2
eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2)
& +akct*deltad*deltat12
- & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
+ & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi+ebr
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
& '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
cosfac=dsqrt(cosfac2)
sinfac2=0.5d0/(1.0d0-costtab(i+1))
sinfac=dsqrt(sinfac2)
- it=itype(i)
+ it=iabs(itype(i))
if (it.eq.10) goto 1
c
C Compute the axes of tghe local cartesian coordinates system; store in
Cc diagnostics - remove later
xx1 = dcos(alph(2))
yy1 = dsin(alph(2))*dcos(omeg(2))
- zz1 = -dsign(1.0,itype(i))*dsin(alph(2))*dsin(omeg(2))
+ zz1 = -dsign(1.0,dfloat(itype(i)))*dsin(alph(2))*dsin(omeg(2))
write(2,'(3f8.1,3f9.3,1x,3f9.3)')
& alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,
& xx1,yy1,zz1
c & dscp1,dscp2,sumene
c sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
escloc = escloc + sumene
-c write (2,*) "i",i," escloc",sumene,escloc
+c write (2,*) "i",i," escloc",sumene,escloc,it,itype(i)
+c & ,zz,xx,yy
c#define DEBUG
#ifdef DEBUG
C
do iii=1,2
dipi(iii,1)=Ub2(iii,i)
dipderi(iii)=Ub2der(iii,i)
- dipi(iii,2)=b1(iii,iti1)
+ dipi(iii,2)=b1(iii,i+1)
dipj(iii,1)=Ub2(iii,j)
dipderj(iii)=Ub2der(iii,j)
- dipj(iii,2)=b1(iii,itj1)
+ dipj(iii,2)=b1(iii,j+1)
enddo
kkk=0
do iii=1,2
C indluded.
IF (wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) THEN
call transpose2(AEA(1,1,1),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,iti),AEAb1(1,1,1))
+ call matvec2(auxmat(1,1),b1(1,i),AEAb1(1,1,1))
call matvec2(auxmat(1,1),Ub2(1,i),AEAb2(1,1,1))
call matvec2(auxmat(1,1),Ub2der(1,i),AEAb2derg(1,2,1,1))
call transpose2(AEAderg(1,1,1),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,iti),AEAb1derg(1,1,1))
+ call matvec2(auxmat(1,1),b1(1,i),AEAb1derg(1,1,1))
call matvec2(auxmat(1,1),Ub2(1,i),AEAb2derg(1,1,1,1))
- call matvec2(AEA(1,1,1),b1(1,itk1),AEAb1(1,2,1))
- call matvec2(AEAderg(1,1,1),b1(1,itk1),AEAb1derg(1,2,1))
+ call matvec2(AEA(1,1,1),b1(1,k+1),AEAb1(1,2,1))
+ call matvec2(AEAderg(1,1,1),b1(1,k+1),AEAb1derg(1,2,1))
call matvec2(AEA(1,1,1),Ub2(1,k+1),AEAb2(1,2,1))
call matvec2(AEAderg(1,1,1),Ub2(1,k+1),AEAb2derg(1,1,2,1))
call matvec2(AEA(1,1,1),Ub2der(1,k+1),AEAb2derg(1,2,2,1))
call transpose2(AEA(1,1,2),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,itj),AEAb1(1,1,2))
+ call matvec2(auxmat(1,1),b1(1,j),AEAb1(1,1,2))
call matvec2(auxmat(1,1),Ub2(1,j),AEAb2(1,1,2))
call matvec2(auxmat(1,1),Ub2der(1,j),AEAb2derg(1,2,1,2))
call transpose2(AEAderg(1,1,2),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,itj),AEAb1derg(1,1,2))
+ call matvec2(auxmat(1,1),b1(1,j),AEAb1derg(1,1,2))
call matvec2(auxmat(1,1),Ub2(1,j),AEAb2derg(1,1,1,2))
- call matvec2(AEA(1,1,2),b1(1,itl1),AEAb1(1,2,2))
- call matvec2(AEAderg(1,1,2),b1(1,itl1),AEAb1derg(1,2,2))
+ call matvec2(AEA(1,1,2),b1(1,l+1),AEAb1(1,2,2))
+ call matvec2(AEAderg(1,1,2),b1(1,l+1),AEAb1derg(1,2,2))
call matvec2(AEA(1,1,2),Ub2(1,l+1),AEAb2(1,2,2))
call matvec2(AEAderg(1,1,2),Ub2(1,l+1),AEAb2derg(1,1,2,2))
call matvec2(AEA(1,1,2),Ub2der(1,l+1),AEAb2derg(1,2,2,2))
do kkk=1,5
do lll=1,3
call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,iti),
+ call matvec2(auxmat(1,1),b1(1,i),
& AEAb1derx(1,lll,kkk,iii,1,1))
call matvec2(auxmat(1,1),Ub2(1,i),
& AEAb2derx(1,lll,kkk,iii,1,1))
- call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,itk1),
+ call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,k+1),
& AEAb1derx(1,lll,kkk,iii,2,1))
call matvec2(AEAderx(1,1,lll,kkk,iii,1),Ub2(1,k+1),
& AEAb2derx(1,lll,kkk,iii,2,1))
call transpose2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,itj),
+ call matvec2(auxmat(1,1),b1(1,j),
& AEAb1derx(1,lll,kkk,iii,1,2))
call matvec2(auxmat(1,1),Ub2(1,j),
& AEAb2derx(1,lll,kkk,iii,1,2))
- call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,itl1),
+ call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,l+1),
& AEAb1derx(1,lll,kkk,iii,2,2))
call matvec2(AEAderx(1,1,lll,kkk,iii,2),Ub2(1,l+1),
& AEAb2derx(1,lll,kkk,iii,2,2))
IF (wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0 .or.
& (wturn6.gt.0.0d0 .and. j.eq.i+4 .and. l.eq.i+3)) THEN
call transpose2(AEA(1,1,1),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,iti),AEAb1(1,1,1))
+ call matvec2(auxmat(1,1),b1(1,i),AEAb1(1,1,1))
call matvec2(auxmat(1,1),Ub2(1,i),AEAb2(1,1,1))
call matvec2(auxmat(1,1),Ub2der(1,i),AEAb2derg(1,2,1,1))
call transpose2(AEAderg(1,1,1),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,iti),AEAb1derg(1,1,1))
+ call matvec2(auxmat(1,1),b1(1,i),AEAb1derg(1,1,1))
call matvec2(auxmat(1,1),Ub2(1,i),AEAb2derg(1,1,1,1))
- call matvec2(AEA(1,1,1),b1(1,itk1),AEAb1(1,2,1))
- call matvec2(AEAderg(1,1,1),b1(1,itk1),AEAb1derg(1,2,1))
+ call matvec2(AEA(1,1,1),b1(1,k+1),AEAb1(1,2,1))
+ call matvec2(AEAderg(1,1,1),b1(1,k+1),AEAb1derg(1,2,1))
call matvec2(AEA(1,1,1),Ub2(1,k+1),AEAb2(1,2,1))
call matvec2(AEAderg(1,1,1),Ub2(1,k+1),AEAb2derg(1,1,2,1))
call matvec2(AEA(1,1,1),Ub2der(1,k+1),AEAb2derg(1,2,2,1))
call transpose2(AEA(1,1,2),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,itj1),AEAb1(1,1,2))
+ call matvec2(auxmat(1,1),b1(1,j+1),AEAb1(1,1,2))
call matvec2(auxmat(1,1),Ub2(1,l),AEAb2(1,1,2))
call matvec2(auxmat(1,1),Ub2der(1,l),AEAb2derg(1,2,1,2))
call transpose2(AEAderg(1,1,2),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,itl),AEAb1(1,1,2))
+ call matvec2(auxmat(1,1),b1(1,l),AEAb1(1,1,2))
call matvec2(auxmat(1,1),Ub2(1,l),AEAb2derg(1,1,1,2))
- call matvec2(AEA(1,1,2),b1(1,itj1),AEAb1(1,2,2))
- call matvec2(AEAderg(1,1,2),b1(1,itj1),AEAb1derg(1,2,2))
+ call matvec2(AEA(1,1,2),b1(1,j+1),AEAb1(1,2,2))
+ call matvec2(AEAderg(1,1,2),b1(1,j+1),AEAb1derg(1,2,2))
call matvec2(AEA(1,1,2),Ub2(1,j),AEAb2(1,2,2))
call matvec2(AEAderg(1,1,2),Ub2(1,j),AEAb2derg(1,1,2,2))
call matvec2(AEA(1,1,2),Ub2der(1,j),AEAb2derg(1,2,2,2))
do kkk=1,5
do lll=1,3
call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,iti),
+ call matvec2(auxmat(1,1),b1(1,i),
& AEAb1derx(1,lll,kkk,iii,1,1))
call matvec2(auxmat(1,1),Ub2(1,i),
& AEAb2derx(1,lll,kkk,iii,1,1))
- call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,itk1),
+ call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,k+1),
& AEAb1derx(1,lll,kkk,iii,2,1))
call matvec2(AEAderx(1,1,lll,kkk,iii,1),Ub2(1,k+1),
& AEAb2derx(1,lll,kkk,iii,2,1))
call transpose2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,itl),
+ call matvec2(auxmat(1,1),b1(1,l),
& AEAb1derx(1,lll,kkk,iii,1,2))
call matvec2(auxmat(1,1),Ub2(1,l),
& AEAb2derx(1,lll,kkk,iii,1,2))
- call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,itj1),
+ call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,j+1),
& AEAb1derx(1,lll,kkk,iii,2,2))
call matvec2(AEAderx(1,1,lll,kkk,iii,2),Ub2(1,j),
& AEAb2derx(1,lll,kkk,iii,2,2))
call matmat2(auxmat(1,1),AEA(1,1,1),pizda(1,1))
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
- eello5_2=scalar2(AEAb1(1,2,1),b1(1,itk))
+ eello5_2=scalar2(AEAb1(1,2,1),b1(1,k))
& -0.5d0*scalar2(vv(1),Ctobr(1,k))
C Explicit gradient in virtual-dihedral angles.
g_corr5_loc(k-1)=g_corr5_loc(k-1)
vv(2)=pizda(2,1)-pizda(1,2)
if (l.eq.j+1) then
g_corr5_loc(l-1)=g_corr5_loc(l-1)
- & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,itk))
+ & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,k))
& -0.5d0*scalar2(vv(1),Ctobr(1,k)))
else
g_corr5_loc(j-1)=g_corr5_loc(j-1)
- & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,itk))
+ & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,k))
& -0.5d0*scalar2(vv(1),Ctobr(1,k)))
endif
C Cartesian gradient
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
derx(lll,kkk,iii)=derx(lll,kkk,iii)
- & +scalar2(AEAb1derx(1,lll,kkk,iii,2,1),b1(1,itk))
+ & +scalar2(AEAb1derx(1,lll,kkk,iii,2,1),b1(1,k))
& -0.5d0*scalar2(vv(1),Ctobr(1,k))
enddo
enddo
call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1))
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
- eello5_4=scalar2(AEAb1(1,2,2),b1(1,itl))
+ eello5_4=scalar2(AEAb1(1,2,2),b1(1,l))
& -0.5d0*scalar2(vv(1),Ctobr(1,l))
C Explicit gradient in virtual-dihedral angles.
g_corr5_loc(l-1)=g_corr5_loc(l-1)
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
g_corr5_loc(k-1)=g_corr5_loc(k-1)
- & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,itl))
+ & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,l))
& -0.5d0*scalar2(vv(1),Ctobr(1,l)))
C Cartesian gradient
do iii=1,2
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
derx(lll,kkk,iii)=derx(lll,kkk,iii)
- & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,itl))
+ & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,l))
& -0.5d0*scalar2(vv(1),Ctobr(1,l))
enddo
enddo
call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1))
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
- eello5_4=scalar2(AEAb1(1,2,2),b1(1,itj))
+ eello5_4=scalar2(AEAb1(1,2,2),b1(1,j))
& -0.5d0*scalar2(vv(1),Ctobr(1,j))
C Explicit gradient in virtual-dihedral angles.
g_corr5_loc(j-1)=g_corr5_loc(j-1)
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
g_corr5_loc(k-1)=g_corr5_loc(k-1)
- & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,itj))
+ & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,j))
& -0.5d0*scalar2(vv(1),Ctobr(1,j)))
C Cartesian gradient
do iii=1,2
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii)
- & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,itj))
+ & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,j))
& -0.5d0*scalar2(vv(1),Ctobr(1,j))
enddo
enddo
vv1(1)=pizda1(1,1)-pizda1(2,2)
vv1(2)=pizda1(1,2)+pizda1(2,1)
s4=0.5d0*scalar2(vv1(1),Dtobr2(1,i))
- vv(1)=AEAb1(1,2,imat)*b1(1,itk)-AEAb1(2,2,imat)*b1(2,itk)
- vv(2)=AEAb1(1,2,imat)*b1(2,itk)+AEAb1(2,2,imat)*b1(1,itk)
+ vv(1)=AEAb1(1,2,imat)*b1(1,k)-AEAb1(2,2,imat)*b1(2,k)
+ vv(2)=AEAb1(1,2,imat)*b1(2,k)+AEAb1(2,2,imat)*b1(1,k)
s5=scalar2(vv(1),Dtobr2(1,i))
cd write (2,*) 's1',s1,' s2',s2,' s3',s3,' s4', s4,' s5',s5
eello6_graph1=-0.5d0*(s1+s2+s3+s4+s5)
call matmat2(AEAderg(1,1,imat),auxmat(1,1),pizda1(1,1))
vv1(1)=pizda1(1,1)-pizda1(2,2)
vv1(2)=pizda1(1,2)+pizda1(2,1)
- vv(1)=AEAb1derg(1,2,imat)*b1(1,itk)-AEAb1derg(2,2,imat)*b1(2,itk)
- vv(2)=AEAb1derg(1,2,imat)*b1(2,itk)+AEAb1derg(2,2,imat)*b1(1,itk)
+ vv(1)=AEAb1derg(1,2,imat)*b1(1,k)-AEAb1derg(2,2,imat)*b1(2,k)
+ vv(2)=AEAb1derg(1,2,imat)*b1(2,k)+AEAb1derg(2,2,imat)*b1(1,k)
if (l.eq.j+1) then
g_corr6_loc(l-1)=g_corr6_loc(l-1)
& +ekont*(-0.5d0*(scalar2(AEAb1derg(1,2,imat),CUgb2(1,i))
vv1(1)=pizda1(1,1)-pizda1(2,2)
vv1(2)=pizda1(1,2)+pizda1(2,1)
s4=0.5d0*scalar2(vv1(1),Dtobr2(1,i))
- vv(1)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(1,itk)
- & -AEAb1derx(2,lll,kkk,iii,2,imat)*b1(2,itk)
- vv(2)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(2,itk)
- & +AEAb1derx(2,lll,kkk,iii,2,imat)*b1(1,itk)
+ vv(1)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(1,k)
+ & -AEAb1derx(2,lll,kkk,iii,2,imat)*b1(2,k)
+ vv(2)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(2,k)
+ & +AEAb1derx(2,lll,kkk,iii,2,imat)*b1(1,k)
s5=scalar2(vv(1),Dtobr2(1,i))
derx(lll,kkk,ind)=derx(lll,kkk,ind)-0.5d0*(s1+s2+s3+s4+s5)
enddo
#ifdef MOMENT
s1=dip(4,jj,i)*dip(4,kk,k)
#endif
- call matvec2(AECA(1,1,1),b1(1,itk1),auxvec(1))
- s2=0.5d0*scalar2(b1(1,itk),auxvec(1))
- call matvec2(AECA(1,1,2),b1(1,itl1),auxvec(1))
- s3=0.5d0*scalar2(b1(1,itj1),auxvec(1))
+ call matvec2(AECA(1,1,1),b1(1,k+1),auxvec(1))
+ s2=0.5d0*scalar2(b1(1,k),auxvec(1))
+ call matvec2(AECA(1,1,2),b1(1,l+1),auxvec(1))
+ s3=0.5d0*scalar2(b1(1,j+1),auxvec(1))
call transpose2(EE(1,1,itk),auxmat(1,1))
call matmat2(auxmat(1,1),AECA(1,1,1),pizda(1,1))
vv(1)=pizda(1,1)+pizda(2,2)
#endif
c eello6_graph3=-s4
C Derivatives in gamma(k-1)
- call matvec2(AECAderg(1,1,2),b1(1,itl1),auxvec(1))
- s3=0.5d0*scalar2(b1(1,itj1),auxvec(1))
+ call matvec2(AECAderg(1,1,2),b1(1,l+1),auxvec(1))
+ s3=0.5d0*scalar2(b1(1,j+1),auxvec(1))
s4=-0.25d0*scalar2(vv(1),Ctobrder(1,k))
g_corr6_loc(k-1)=g_corr6_loc(k-1)-ekont*(s3+s4)
C Derivatives in gamma(l-1)
- call matvec2(AECAderg(1,1,1),b1(1,itk1),auxvec(1))
- s2=0.5d0*scalar2(b1(1,itk),auxvec(1))
+ call matvec2(AECAderg(1,1,1),b1(1,k+1),auxvec(1))
+ s2=0.5d0*scalar2(b1(1,k),auxvec(1))
call matmat2(auxmat(1,1),AECAderg(1,1,1),pizda(1,1))
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
s1=dip(4,jj,i)*dipderx(lll,kkk,4,kk,k)
endif
#endif
- call matvec2(AECAderx(1,1,lll,kkk,iii,1),b1(1,itk1),
+ call matvec2(AECAderx(1,1,lll,kkk,iii,1),b1(1,k+1),
& auxvec(1))
- s2=0.5d0*scalar2(b1(1,itk),auxvec(1))
- call matvec2(AECAderx(1,1,lll,kkk,iii,2),b1(1,itl1),
+ s2=0.5d0*scalar2(b1(1,k),auxvec(1))
+ call matvec2(AECAderx(1,1,lll,kkk,iii,2),b1(1,l+1),
& auxvec(1))
- s3=0.5d0*scalar2(b1(1,itj1),auxvec(1))
+ s3=0.5d0*scalar2(b1(1,j+1),auxvec(1))
call matmat2(auxmat(1,1),AECAderx(1,1,lll,kkk,iii,1),
& pizda(1,1))
vv(1)=pizda(1,1)+pizda(2,2)
call matvec2(AECA(1,1,imat),Ub2(1,k),auxvec(1))
s2=0.5d0*scalar2(Ub2(1,i),auxvec(1))
if (j.eq.l+1) then
- call matvec2(ADtEA1(1,1,3-imat),b1(1,itj1),auxvec1(1))
- s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1))
+ call matvec2(ADtEA1(1,1,3-imat),b1(1,j+1),auxvec1(1))
+ s3=-0.5d0*scalar2(b1(1,j),auxvec1(1))
else
- call matvec2(ADtEA1(1,1,3-imat),b1(1,itl1),auxvec1(1))
- s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1))
+ call matvec2(ADtEA1(1,1,3-imat),b1(1,l+1),auxvec1(1))
+ s3=-0.5d0*scalar2(b1(1,l),auxvec1(1))
endif
call transpose2(EUg(1,1,k),auxmat(1,1))
call matmat2(AECA(1,1,imat),auxmat(1,1),pizda(1,1))
#endif
s2=0.5d0*scalar2(Ub2der(1,i),auxvec(1))
if (j.eq.l+1) then
- call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,itj1),auxvec1(1))
- s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1))
+ call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,j+1),auxvec1(1))
+ s3=-0.5d0*scalar2(b1(1,j),auxvec1(1))
else
- call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,itl1),auxvec1(1))
- s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1))
+ call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,l+1),auxvec1(1))
+ s3=-0.5d0*scalar2(b1(1,l),auxvec1(1))
endif
s4=0.25d0*scalar2(vv(1),Dtobr2der(1,i))
if (wturn6.gt.0.0d0 .and. k.eq.l+4 .and. i.eq.j+2) then
call matvec2(AECA(1,1,imat),Ub2der(1,k),auxvec1(1))
s2=0.5d0*scalar2(Ub2(1,i),auxvec1(1))
if (j.eq.l+1) then
- call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,itj1),auxvec1(1))
- s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1))
+ call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,j+1),auxvec1(1))
+ s3=-0.5d0*scalar2(b1(1,j),auxvec1(1))
else
- call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,itl1),auxvec1(1))
- s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1))
+ call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,l+1),auxvec1(1))
+ s3=-0.5d0*scalar2(b1(1,l),auxvec1(1))
endif
call transpose2(EUgder(1,1,k),auxmat1(1,1))
call matmat2(AECA(1,1,imat),auxmat1(1,1),pizda(1,1))
s2=0.5d0*scalar2(Ub2(1,i),auxvec(1))
if (j.eq.l+1) then
call matvec2(ADtEA1derx(1,1,lll,kkk,iii,3-imat),
- & b1(1,itj1),auxvec(1))
- s3=-0.5d0*scalar2(b1(1,itj),auxvec(1))
+ & b1(1,j+1),auxvec(1))
+ s3=-0.5d0*scalar2(b1(1,j),auxvec(1))
else
call matvec2(ADtEA1derx(1,1,lll,kkk,iii,3-imat),
- & b1(1,itl1),auxvec(1))
- s3=-0.5d0*scalar2(b1(1,itl),auxvec(1))
+ & b1(1,l+1),auxvec(1))
+ s3=-0.5d0*scalar2(b1(1,l),auxvec(1))
endif
call matmat2(AECAderx(1,1,lll,kkk,iii,imat),auxmat(1,1),
& pizda(1,1))
#ifdef MOMENT
call transpose2(AEA(1,1,1),auxmat(1,1))
call matmat2(EUg(1,1,i+1),auxmat(1,1),auxmat(1,1))
- ss1=scalar2(Ub2(1,i+2),b1(1,itl))
+ ss1=scalar2(Ub2(1,i+2),b1(1,l))
s1 = (auxmat(1,1)+auxmat(2,2))*ss1
#endif
- call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1(1))
+ call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1(1))
call matvec2(AEA(1,1,1),vtemp1(1),vtemp1(1))
- s2 = scalar2(b1(1,itk),vtemp1(1))
+ s2 = scalar2(b1(1,k),vtemp1(1))
#ifdef MOMENT
call transpose2(AEA(1,1,2),atemp(1,1))
call matmat2(atemp(1,1),EUg(1,1,i+4),atemp(1,1))
call matmat2(achuj_temp(1,1),EUg(1,1,i+2),gtemp(1,1))
call matmat2(gtemp(1,1),EUg(1,1,i+3),gtemp(1,1))
call matvec2(a_chuj(1,1,jj,i),Ub2(1,i+4),vtemp4(1))
- ss13 = scalar2(b1(1,itk),vtemp4(1))
+ ss13 = scalar2(b1(1,k),vtemp4(1))
s13 = (gtemp(1,1)+gtemp(2,2))*ss13
#endif
c write (2,*) 's1,s2,s8,s12,s13',s1,s2,s8,s12,s13
#ifdef MOMENT
call transpose2(AEA(1,1,1),auxmatd(1,1))
call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1))
- ss1d=scalar2(Ub2der(1,i+2),b1(1,itl))
+ ss1d=scalar2(Ub2der(1,i+2),b1(1,l))
s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1d
#endif
- call matvec2(EUgder(1,1,i+2),b1(1,itl),vtemp1d(1))
+ call matvec2(EUgder(1,1,i+2),b1(1,l),vtemp1d(1))
call matvec2(AEA(1,1,1),vtemp1d(1),vtemp1d(1))
- s2d = scalar2(b1(1,itk),vtemp1d(1))
+ s2d = scalar2(b1(1,k),vtemp1d(1))
#ifdef MOMENT
call matvec2(Ug2der(1,1,i+2),dd(1,1,itk1),vtemp2d(1))
s8d = -(atemp(1,1)+atemp(2,2))*scalar2(cc(1,1,itl),vtemp2d(1))
call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1))
s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1
#endif
- call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1d(1))
+ call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1d(1))
call matvec2(AEAderg(1,1,1),vtemp1d(1),vtemp1d(1))
- s2d = scalar2(b1(1,itk),vtemp1d(1))
+ s2d = scalar2(b1(1,k),vtemp1d(1))
#ifdef MOMENT
call transpose2(AEA(1,1,2),atempd(1,1))
call matmat2(atempd(1,1),EUgder(1,1,i+4),atempd(1,1))
s12d = scalar2(Ub2(1,i+2),vtemp3d(1))
#ifdef MOMENT
call matvec2(a_chuj(1,1,jj,i),Ub2der(1,i+4),vtemp4d(1))
- ss13d = scalar2(b1(1,itk),vtemp4d(1))
+ ss13d = scalar2(b1(1,k),vtemp4d(1))
s13d = (gtemp(1,1)+gtemp(2,2))*ss13d
#endif
c s1d=0.0d0
call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1))
s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1
#endif
- call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1(1))
+ call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1(1))
call matvec2(AEAderx(1,1,lll,kkk,iii,1),vtemp1(1),
& vtemp1d(1))
- s2d = scalar2(b1(1,itk),vtemp1d(1))
+ s2d = scalar2(b1(1,k),vtemp1d(1))
#ifdef MOMENT
call transpose2(AEAderx(1,1,lll,kkk,iii,2),atempd(1,1))
call matmat2(atempd(1,1),EUg(1,1,i+4),atempd(1,1))
derx_turn(lll,kkk,2) = derx_turn(lll,kkk,2)-0.5d0*s13d
call matvec2(a_chuj_der(1,1,lll,kkk,jj,i),Ub2(1,i+4),
& vtemp4d(1))
- ss13d = scalar2(b1(1,itk),vtemp4d(1))
+ ss13d = scalar2(b1(1,k),vtemp4d(1))
s13d = (gtemp(1,1)+gtemp(2,2))*ss13d
derx_turn(lll,kkk,1) = derx_turn(lll,kkk,1)-0.5d0*s13d
enddo