X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-M%2Fenergy_p_new.F;h=cac8518ddb97a2ee4fd112196b6958f1954e7c61;hb=6f877b0a7b5b4f292dd52b2b447ececb7e469044;hp=b266ab4a21655d9d23ccb51e8d8f164c11cd2342;hpb=a12fec1b1bb1405b29ce7db86524648c70f5d635;p=unres.git diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index b266ab4..cac8518 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -1916,10 +1916,15 @@ C 12/26/95 - for the evaluation of multi-body H-bonding interactions ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg) ees=ees+eesij evdw1=evdw1+evdwij -cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') -cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, -cd & 1.0D0/dsqrt(rrmij),evdwij,eesij, -cd & xmedi,ymedi,zmedi,xj,yj,zj +c write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') +c &'evdw1',i,j,evdwij +c &,iteli,itelj,aaa,evdw1 + +c write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij +c write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') +c & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, +c & 1.0D0/dsqrt(rrmij),evdwij,eesij, +c & xmedi,ymedi,zmedi,xj,yj,zj C C Calculate contributions to the Cartesian gradient. C @@ -2265,8 +2270,10 @@ C Check the loc-el terms by numerical integration 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) -cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij -cd write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) +c write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij +c write (iout,'(a6,2i5,0pf7.3)') +c & '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 if (calc_grad) then @@ -2835,9 +2842,9 @@ C Uncomment following three lines for Ca-p interactions evdw2_14=evdw2_14+e1+e2 endif evdwij=e1+e2 - write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)') - & 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli), - & bad(itypj,iteli) +c write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)') +c & 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli), +c & bad(itypj,iteli) evdw2=evdw2+evdwij if (calc_grad) then C @@ -3172,8 +3179,8 @@ C & delthe0,sig0inv,sigtc,sigsqtc,delthec,it double precision y(2),z(2) delta=0.02d0*pi - time11=dexp(-2*time) - time12=1.0d0 +c time11=dexp(-2*time) +c time12=1.0d0 etheta=0.0D0 c write (iout,*) "nres",nres c write (*,'(a,i2)') 'EBEND ICG=',icg @@ -3199,8 +3206,8 @@ C Zero the energy function and its derivative at 0 or pi. if (i.gt.3 .and. itype(i-2).ne.ntyp1) then #ifdef OSF phii=phi(i) - icrc=0 - call proc_proc(phii,icrc) +c icrc=0 +c call proc_proc(phii,icrc) if (icrc.eq.1) phii=150.0 #else phii=phi(i) @@ -3214,8 +3221,8 @@ C Zero the energy function and its derivative at 0 or pi. if (i.lt.nres .and. itype(i).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) - icrc=0 - call proc_proc(phii1,icrc) +c icrc=0 +c call proc_proc(phii1,icrc) if (icrc.eq.1) phii1=150.0 phii1=pinorm(phii1) z(1)=cos(phii1) @@ -3283,7 +3290,7 @@ c & rad2deg*phii,rad2deg*phii1,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) - 1215 continue +c 1215 continue enddo C Ufff.... We've done all this!!! return @@ -3426,7 +3433,9 @@ C etheta=0.0D0 c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) do i=ithet_start,ithet_end - if (itype(i-1).eq.ntyp1) cycle +c if (itype(i-1).eq.ntyp1) cycle + if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. + &(itype(i).eq.ntyp1)) cycle if (iabs(itype(i+1)).eq.20) iblock=2 if (iabs(itype(i+1)).ne.20) iblock=1 dethetai=0.0d0 @@ -3438,7 +3447,7 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3 .and. itype(i-2).ne.ntyp1) then + if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -3452,13 +3461,14 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) enddo else phii=0.0d0 - ityp1=nthetyp+1 +c ityp1=nthetyp+1 do k=1,nsingle + ityp1=ithetyp((itype(i-2))) cosph1(k)=0.0d0 sinph1(k)=0.0d0 enddo endif - if (i.lt.nres .and. itype(i).ne.ntyp1) then + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -3473,7 +3483,8 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) enddo else phii1=0.0d0 - ityp3=nthetyp+1 +c ityp3=nthetyp+1 + ityp3=ithetyp((itype(i))) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 @@ -3590,7 +3601,8 @@ c call flush(iout) 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 +c gloc(nphi+i-2,icg)=wang*dethetai + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai enddo return end @@ -3941,7 +3953,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)*dsign(1.0,dfloat(itype(i))) + 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) @@ -3981,7 +3993,7 @@ c write (2,*) "xx",xx," yy",yy," zz",zz 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.0d0,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 @@ -4024,6 +4036,8 @@ c & dscp1,dscp2,sumene c sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1)) escloc = escloc + sumene c write (2,*) "escloc",escloc +c write (2,*) "i",i," escloc",sumene,escloc,it,itype(i), +c & zz,xx,yy if (.not. calc_grad) goto 1 #ifdef DEBUG C @@ -4153,9 +4167,9 @@ c & (dC_norm(j,i-1),j=1,3)," vbld_inv",vbld_inv(i+1),vbld_inv(i) dZZ_Ci(k)=0.0d0 do j=1,3 dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1) - & *dsign(1.0,dfloat(itype(i)))*dC_norm(j,i+nres) + & *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.0,dfloat(itype(i)))*dC_norm(j,i+nres) + & *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)) @@ -4726,9 +4740,9 @@ c------------------------------------------------------------------------------ integer dimen1,dimen2,atom,indx double precision buffer(dimen1,dimen2) double precision zapas - common /contacts_hb/ zapas(3,20,maxres,7), - & facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres), - & num_cont_hb(maxres),jcont_hb(20,maxres) + common /contacts_hb/ zapas(3,ntyp,maxres,7), + & facont_hb(ntyp,maxres),ees0p(ntyp,maxres),ees0m(ntyp,maxres), + & num_cont_hb(maxres),jcont_hb(ntyp,maxres) num_kont=num_cont_hb(atom) do i=1,num_kont do k=1,7 @@ -4751,9 +4765,10 @@ c------------------------------------------------------------------------------ integer dimen1,dimen2,atom,indx double precision buffer(dimen1,dimen2) double precision zapas - common /contacts_hb/ zapas(3,20,maxres,7), - & facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres), - & num_cont_hb(maxres),jcont_hb(20,maxres) + common /contacts_hb/ zapas(3,ntyp,maxres,7), + & facont_hb(ntyp,maxres),ees0p(ntyp,maxres), + & ees0m(ntyp,maxres), + & num_cont_hb(maxres),jcont_hb(ntyp,maxres) num_kont=buffer(1,indx+26) num_kont_old=num_cont_hb(atom) num_cont_hb(atom)=num_kont+num_kont_old @@ -6492,7 +6507,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