X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=97986d218d5a50790533392449a571059a087efb;hb=153cbfdee6e22c0c33798211afef1ec4cf01c117;hp=38c530654bb79276f87b5381a323411ef336a41e;hpb=1fd1d1d47d5af85ec47d73bede0f98e38ca1f32a;p=unres.git 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 38c5306..97986d2 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -2938,7 +2938,9 @@ cd & 1.0D0/dsqrt(rrmij),evdwij,eesij, 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 @@ -3956,8 +3958,9 @@ C Uncomment following three lines for Ca-p interactions endif evdwij=e1+e2 evdw2=evdw2+evdwij - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') - & 'evdw2',i,j,evdwij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)') + & 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli), + & bad(itypj,iteli) C C Calculate contributions to the gradient in the virtual-bond and SC vectors. C @@ -4545,11 +4548,13 @@ C etheta=0.0D0 do i=ithet_start,ithet_end if (itype(i-1).eq.ntyp1) cycle + if (iabs(itype(i+1)).eq.20) iblock=2 + if (iabs(itype(i+1)).ne.20) iblock=1 dethetai=0.0d0 dephii=0.0d0 dephii1=0.0d0 theti2=0.5d0*theta(i) - ityp2=ithetyp(iabs(itype(i-1))) + ityp2=ithetyp((itype(i-1))) do k=1,nntheterm coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) @@ -4561,7 +4566,8 @@ C #else phii=phi(i) #endif - ityp1=ithetyp(iabs(itype(i-2))) + ityp1=ithetyp((itype(i-2))) +C propagation of chirality for glycine type do k=1,nsingle cosph1(k)=dcos(k*phii) sinph1(k)=dsin(k*phii) @@ -4582,7 +4588,7 @@ C #else phii1=phi(i+1) #endif - ityp3=ithetyp(iabs(itype(i))) + ityp3=ithetyp((itype(i))) do k=1,nsingle cosph2(k)=dcos(k*phii1) sinph2(k)=dsin(k*phii1) @@ -4595,7 +4601,7 @@ C sinph2(k)=0.0d0 enddo endif - ethetai=aa0thet(ityp1,ityp2,ityp3) + ethetai=aa0thet(ityp1,ityp2,ityp3,iblock) do k=1,ndouble do l=1,k-1 ccl=cosph1(l)*cosph2(k-l) @@ -4617,11 +4623,12 @@ C enddo endif do k=1,ntheterm - ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3)*sinkt(k) - dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3) + ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3,iblock)*sinkt(k) + dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3,iblock) & *coskt(k) if (lprn) - & write (iout,*) "k",k," aathet",aathet(k,ityp1,ityp2,ityp3), + & write (iout,*) "k",k," + & aathet",aathet(k,ityp1,ityp2,ityp3,iblock), & " ethetai",ethetai enddo if (lprn) then @@ -4640,24 +4647,24 @@ C endif do m=1,ntheterm2 do k=1,nsingle - aux=bbthet(k,m,ityp1,ityp2,ityp3)*cosph1(k) - & +ccthet(k,m,ityp1,ityp2,ityp3)*sinph1(k) - & +ddthet(k,m,ityp1,ityp2,ityp3)*cosph2(k) - & +eethet(k,m,ityp1,ityp2,ityp3)*sinph2(k) + aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k) + & +ccthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k) + & +ddthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k) + & +eethet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k) ethetai=ethetai+sinkt(m)*aux dethetai=dethetai+0.5d0*m*aux*coskt(m) dephii=dephii+k*sinkt(m)*( - & ccthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)- - & bbthet(k,m,ityp1,ityp2,ityp3)*sinph1(k)) + & ccthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)- + & bbthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k)) dephii1=dephii1+k*sinkt(m)*( - & eethet(k,m,ityp1,ityp2,ityp3)*cosph2(k)- - & ddthet(k,m,ityp1,ityp2,ityp3)*sinph2(k)) + & eethet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)- + & ddthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k)) if (lprn) & write (iout,*) "m",m," k",k," bbthet", - & bbthet(k,m,ityp1,ityp2,ityp3)," ccthet", - & ccthet(k,m,ityp1,ityp2,ityp3)," ddthet", - & ddthet(k,m,ityp1,ityp2,ityp3)," eethet", - & eethet(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai + & bbthet(k,m,ityp1,ityp2,ityp3,iblock)," ccthet", + & ccthet(k,m,ityp1,ityp2,ityp3,iblock)," ddthet", + & ddthet(k,m,ityp1,ityp2,ityp3,iblock)," eethet", + & eethet(k,m,ityp1,ityp2,ityp3,iblock)," ethetai",ethetai enddo enddo if (lprn) @@ -4665,28 +4672,29 @@ C do m=1,ntheterm3 do k=2,ndouble do l=1,k-1 - aux=ffthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+ - & ffthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+ - & ggthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+ - & ggthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l) + aux=ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+ + & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l)+ + & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+ + & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l) ethetai=ethetai+sinkt(m)*aux dethetai=dethetai+0.5d0*m*coskt(m)*aux dephii=dephii+l*sinkt(m)*( - & -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)- - & ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+ - & ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+ - & ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)) + & -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)- + & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+ + & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+ + & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l)) dephii1=dephii1+(k-l)*sinkt(m)*( - & -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+ - & ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+ - & ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)- - & ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)) + & -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+ + & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+ + & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)- + & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l)) if (lprn) then write (iout,*) "m",m," k",k," l",l," ffthet", - & ffthet(l,k,m,ityp1,ityp2,ityp3), - & ffthet(k,l,m,ityp1,ityp2,ityp3)," ggthet", - & ggthet(l,k,m,ityp1,ityp2,ityp3), - & ggthet(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai + & ffthet(l,k,m,ityp1,ityp2,ityp3,iblock), + & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)," ggthet", + & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock), + & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock), + & " ethetai",ethetai write (iout,*) cosph1ph2(l,k)*sinkt(m), & cosph1ph2(k,l)*sinkt(m), & sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m) @@ -4695,9 +4703,12 @@ C enddo enddo 10 continue - if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') +c lprn1=.true. + if (lprn1) + & write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') & i,theta(i)*rad2deg,phii*rad2deg, & phii1*rad2deg,ethetai +c lprn1=.false. 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 @@ -5058,7 +5069,7 @@ C & dc_norm(3,i+nres) y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac enddo do j = 1,3 - z_prime(j) = -uz(j,i-1) + z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i))) enddo c write (2,*) "i",i c write (2,*) "x_prime",(x_prime(j),j=1,3) @@ -5080,7 +5091,7 @@ c do j = 1,3 xx = xx + x_prime(j)*dc_norm(j,i+nres) yy = yy + y_prime(j)*dc_norm(j,i+nres) - zz = zz + dsign(1.0,itype(i))*z_prime(j)*dc_norm(j,i+nres) + zz = zz + z_prime(j)*dc_norm(j,i+nres) enddo xxtab(i)=xx @@ -5141,6 +5152,7 @@ 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#define DEBUG #ifdef DEBUG C C This section to check the numerical derivatives of the energy of ith side @@ -5184,6 +5196,7 @@ C End of diagnostics section. C C Compute the gradient of esc C +c zz=zz*dsign(1.0,dfloat(itype(i))) pom_s1=(1.0d0+x(63))/(0.1d0 + dscp1)**2 pom_s16=6*(1.0d0+x(64))/(0.1d0 + dscp1**6)**2 pom_s2=(1.0d0+x(65))/(0.1d0 + dscp2)**2 @@ -5208,7 +5221,7 @@ C & +(sumene2x+sumene4x*cost2tab(i+1))*(s2+s2_6) & +(pom1+pom2)*pom_dx #ifdef DEBUG - write(2,*), "de_dxx = ", de_dxx,de_dxx_num + write(2,*), "de_dxx = ", de_dxx,de_dxx_num,itype(i) #endif C sumene1y=x(3) + 2*x(6)*yy + x(9)*xx + x(10)*zz @@ -5223,7 +5236,7 @@ C & +(sumene2y+sumene4y*cost2tab(i+1))*(s2+s2_6) & +(pom1-pom2)*pom_dy #ifdef DEBUG - write(2,*), "de_dyy = ", de_dyy,de_dyy_num + write(2,*), "de_dyy = ", de_dyy,de_dyy_num,itype(i) #endif C de_dzz =(x(24) +2*x(27)*zz +x(28)*xx +x(30)*yy @@ -5235,15 +5248,16 @@ C & +x(60)*xx*yy)*cost2tab(i+1)*(s2+s2_6) & + ( x(14) + 2*x(17)*zz+ x(18)*xx + x(20)*yy)*(s2+s2_6) #ifdef DEBUG - write(2,*), "de_dzz = ", de_dzz,de_dzz_num + write(2,*), "de_dzz = ", de_dzz,de_dzz_num,itype(i) #endif C de_dt = 0.5d0*sumene3*cost2tab(i+1)*(s1+s1_6) & -0.5d0*sumene4*sint2tab(i+1)*(s2+s2_6) & +pom1*pom_dt1+pom2*pom_dt2 #ifdef DEBUG - write(2,*), "de_dt = ", de_dt,de_dt_num + write(2,*), "de_dt = ", de_dt,de_dt_num,itype(i) #endif +c#undef DEBUG c C cossc=scalar(dc_norm(1,i),dc_norm(1,i+nres)) @@ -5268,13 +5282,16 @@ c & (dC_norm(j,i-1),j=1,3)," vbld_inv",vbld_inv(i+1),vbld_inv(i) dZZ_Ci1(k)=0.0d0 dZZ_Ci(k)=0.0d0 do j=1,3 - dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres) - dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres) + dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1) + & *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres) + dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1) + & *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres) enddo dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres)) dYY_XYZ(k)=vbld_inv(i+nres)*(y_prime(k)-yy*dC_norm(k,i+nres)) - dZZ_XYZ(k)=vbld_inv(i+nres)*(z_prime(k)-zz*dC_norm(k,i+nres)) + dZZ_XYZ(k)=vbld_inv(i+nres)* + & (z_prime(k)-zz*dC_norm(k,i+nres)) c dt_dCi(k) = -dt_dCi(k)/sinttab(i+1) dt_dCi1(k)= -dt_dCi1(k)/sinttab(i+1) @@ -5656,6 +5673,7 @@ C Set lprn=.true. for debugging lprn=.false. c lprn=.true. etors_d=0.0D0 +c write(iout,*) "a tu??" do i=iphid_start,iphid_end if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 & .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle @@ -5737,6 +5755,7 @@ c lprn=.true. c write (iout,*) "EBACK_SC_COR",itau_start,itau_end esccor=0.0D0 do i=itau_start,itau_end + if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle esccor_ii=0.0D0 isccori=isccortyp(itype(i-2)) isccori1=isccortyp(itype(i-1)) @@ -5769,7 +5788,7 @@ c 3 = SC...Ca...Ca...SCi esccor=esccor+v1ij*cosphi+v2ij*sinphi gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) enddo -c write (iout,*) "EBACK_SC_COR",i,esccor,intertyp +c write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci if (lprn) & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') @@ -7983,7 +8002,7 @@ c---------------------------------------------------------------------------- include 'COMMON.GEO' logical swap double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2), - & auxvec1(2),auxvec2(1),auxmat1(2,2) + & auxvec1(2),auxvec2(2),auxmat1(2,2) logical lprn common /kutas/ lprn CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC