From: Adam Sieradzan Date: Thu, 14 Nov 2013 22:25:51 +0000 (+0100) Subject: Dzialajacy gradient dla reszt 13 i 7 X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;h=02826f782138c5757ee627b1785503cc313a1309;p=unres.git Dzialajacy gradient dla reszt 13 i 7 --- diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index c13d341..2682840 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -2252,7 +2252,7 @@ C 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 @@ -2270,23 +2270,23 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then 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, @@ -2307,8 +2307,9 @@ c b2(2,iti)=bnew2(1,2,iti)*sin(alpha(i))*sin(beta(i)) 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 @@ -2447,7 +2448,7 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then 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) @@ -2895,7 +2896,7 @@ 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=7,7 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -2909,8 +2910,8 @@ c do i=iatel_s,iatel_e 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=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 @@ -3181,13 +3182,14 @@ C 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 @@ -3354,37 +3356,47 @@ cgrad endif 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