else
iti1=ntortyp+1
endif
- write(iout,*),i
+c 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(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,*) '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
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
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
+c 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
+c 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
& +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)
+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)
& +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)+
& 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)
+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))
& '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