C Compute the virtual-bond-torsional-angle dependent quantities needed
C to calculate the el-loc multibody terms of various order.
C
+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
+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)*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)*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,
+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)
+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
+#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
iti1=ntortyp
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)
+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)
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
if (zmedi.lt.0) zmedi=zmedi+boxzsize
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)
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
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)
+c write(iout,*) 'mumu=', mu(k,i),mu(l,j),i,j,k,l
+#ifdef NEWCORR
+ 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
cd write (iout,*) 'EELEC: i',i,' j',j
& +a33*muij(4)
c write (iout,*) 'i',i,' j',j,itype(i),itype(j),
c & ' eel_loc_ij',eel_loc_ij
+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)
+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
+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
+c Derivative over j residue
+ geel_loc_ji=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)+
+ & geel_loc_ji*wel_loc
+ geel_loc_ji=
+ & +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
+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
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)
c write(iout,*)'chujOWO', auxvec(1),b1(1,iti2)
if (energy_dec) write (iout,'(a6,2i5,0pf7.3,3f7.3)')
& 'eturn4',i,j,-(s1+s2+s3),s1,s2,s3
cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
cd & ' eello_turn4_num',8*eello_turn4_num
+#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)
+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))
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)
+ gloc(nphi+i-2,icg)=wang*dethetai+gloc(nphi+i-2,icg)
enddo
return
end
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
projects / unres.git / blobdiff