X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-M%2Fenergy_p_new.F;h=7217542b604fb7236557707fe99542d61f5cc2f6;hb=1e3c2bda475e5d4031b8dab868bb8bf379a9e388;hp=00ed59ce122e5019a7dec24acd5c889f6a9094e2;hpb=dbb0c0d7829893e0653c15ad74f2ae8a98bf3ea5;p=unres.git diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index 00ed59c..7217542 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -44,11 +44,13 @@ C Gay-Berne potential (shifted LJ, angular dependence). goto 106 C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence). 105 call egbv(evdw,evdw_t) +C write(iout,*) 'po elektostatyce' C C Calculate electrostatic (H-bonding) energy of the main chain. C - 106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) -C + 106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) +C write(iout,*) 'po eelec' + C Calculate excluded-volume interaction energy between peptide groups C and side chains. C @@ -56,8 +58,9 @@ C c c Calculate the bond-stretching energy c + call ebond(estr) -c write (iout,*) "estr",estr +C write (iout,*) "estr",estr C C Calculate the disulfide-bridge and other energy and the contributions C from other distance constraints. @@ -68,12 +71,12 @@ C C Calculate the virtual-bond-angle energy. C call ebend(ebe) -cd print *,'Bend energy finished.' +C print *,'Bend energy finished.' C C Calculate the SC local energy. C call esc(escloc) -cd print *,'SCLOC energy finished.' +C print *,'SCLOC energy finished.' C C Calculate the virtual-bond torsional energy. C @@ -797,6 +800,14 @@ c if (icall.gt.0) lprn=.true. xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) +C returning the ith atom to box + xi=mod(xi,boxxsize) + if (xi.lt.0) xi=xi+boxxsize + yi=mod(yi,boxysize) + if (yi.lt.0) yi=yi+boxysize + zi=mod(zi,boxzsize) + if (zi.lt.0) zi=zi+boxzsize + dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -830,15 +841,59 @@ c chip12=0.0D0 c alf1=0.0D0 c alf2=0.0D0 c alf12=0.0D0 - xj=c(1,nres+j)-xi - yj=c(2,nres+j)-yi - zj=c(3,nres+j)-zi + xj=c(1,nres+j) + yj=c(2,nres+j) + zj=c(3,nres+j) +C returning jth atom to box + xj=mod(xj,boxxsize) + if (xj.lt.0) xj=xj+boxxsize + yj=mod(yj,boxysize) + if (yj.lt.0) yj=yj+boxysize + zj=mod(zj,boxzsize) + if (zj.lt.0) zj=zj+boxzsize +C checking the distance + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 +C finding the closest + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif + dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) c write (iout,*) i,j,xj,yj,zj rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) + sss=sscale((1.0d0/rij)/sigma(itypi,itypj)) + sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj)) + if (sss.le.0.0) cycle C Calculate angle-dependent terms of energy and contributions to their C derivatives. call sc_angular @@ -861,9 +916,9 @@ c--------------------------------------------------------------- eps3der=evdwij*eps2rt evdwij=evdwij*eps2rt*eps3rt if (bb(itypi,itypj).gt.0) then - evdw=evdw+evdwij + evdw=evdw+evdwij*sss else - evdw_t=evdw_t+evdwij + evdw_t=evdw_t+evdwij*sss endif ij=icant(itypi,itypj) aux=eps1*eps2rt**2*eps3rt**2 @@ -892,6 +947,7 @@ C Calculate gradient components. fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder fac=rij*fac + fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac @@ -1845,14 +1901,22 @@ cd enddo do i=1,nres num_cont_hb(i)=0 enddo -cd print '(a)','Enter EELEC' -cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e +C print '(a)','Enter EELEC' +C write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e do i=1,nres gel_loc_loc(i)=0.0d0 gcorr_loc(i)=0.0d0 enddo do i=iatel_s,iatel_e - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle + if (i.eq.1) then + if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1 + & .or. itype(i+2).eq.ntyp1) cycle + else + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 + & .or. itype(i+2).eq.ntyp1 + & .or. itype(i-1).eq.ntyp1 + &) cycle + endif if (itel(i).eq.0) goto 1215 dxi=dc(1,i) dyi=dc(2,i) @@ -1863,10 +1927,27 @@ cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi + xmedi=mod(xmedi,boxxsize) + if (xmedi.lt.0) xmedi=xmedi+boxxsize + ymedi=mod(ymedi,boxysize) + if (ymedi.lt.0) ymedi=ymedi+boxysize + zmedi=mod(zmedi,boxzsize) + if (zmedi.lt.0) zmedi=zmedi+boxzsize num_conti=0 -c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) +C write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) do j=ielstart(i),ielend(i) - if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle + if (j.eq.1) then + if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 + & .or.itype(j+2).eq.ntyp1 + &) cycle + else + if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 + & .or.itype(j+2).eq.ntyp1 + & .or.itype(j-1).eq.ntyp1 + &) cycle + endif +C +C) cycle if (itel(j).eq.0) goto 1216 ind=ind+1 iteli=itel(i) @@ -1888,10 +1969,50 @@ C End diagnostics dx_normj=dc_norm(1,j) dy_normj=dc_norm(2,j) dz_normj=dc_norm(3,j) - xj=c(1,j)+0.5D0*dxj-xmedi - yj=c(2,j)+0.5D0*dyj-ymedi - zj=c(3,j)+0.5D0*dzj-zmedi + xj=c(1,j)+0.5D0*dxj + yj=c(2,j)+0.5D0*dyj + zj=c(3,j)+0.5D0*dzj + xj=mod(xj,boxxsize) + if (xj.lt.0) xj=xj+boxxsize + yj=mod(yj,boxysize) + if (yj.lt.0) yj=yj+boxysize + zj=mod(zj,boxzsize) + if (zj.lt.0) zj=zj+boxzsize + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + isubchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + isubchap=1 + endif + enddo + enddo + enddo + if (isubchap.eq.1) then + xj=xj_temp-xmedi + yj=yj_temp-ymedi + zj=zj_temp-zmedi + else + xj=xj_safe-xmedi + yj=yj_safe-ymedi + zj=zj_safe-zmedi + endif + rij=xj*xj+yj*yj+zj*zj + sss=sscale(sqrt(rij)) + sssgrad=sscagrad(sqrt(rij)) rrmij=1.0D0/rij rij=dsqrt(rij) rmij=1.0D0/rij @@ -1915,21 +2036,21 @@ c write (iout,*) "i",i,iteli," j",j,itelj," eesij",eesij 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 + evdw1=evdw1+evdwij*sss c write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') c &'evdw1',i,j,evdwij c &,iteli,itelj,aaa,evdw1 - write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij -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)') '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 #ifdef SPLITELE - facvdw=-6*rrmij*(ev1+evdwij) + facvdw=-6*rrmij*(ev1+evdwij)*sss facel=-3*rrmij*(el1+eesij) fac1=fac erij(1)=xj*rmij @@ -1955,9 +2076,18 @@ C gelc(l,k)=gelc(l,k)+ggg(l) enddo enddo - ggg(1)=facvdw*xj - ggg(2)=facvdw*yj - ggg(3)=facvdw*zj +C ggg(1)=facvdw*xj +C ggg(2)=facvdw*yj +C ggg(3)=facvdw*zj + if (sss.gt.0.0) then + ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj + ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj + else + ggg(1)=0.0 + ggg(2)=0.0 + ggg(3)=0.0 + endif do k=1,3 ghalf=0.5D0*ggg(k) gvdwpp(k,i)=gvdwpp(k,i)+ghalf @@ -1972,7 +2102,7 @@ C enddo enddo #else - facvdw=ev1+evdwij + facvdw=(ev1+evdwij)*sss facel=el1+eesij fac1=fac fac=-3*rrmij*(facvdw+facvdw+facel) @@ -2270,8 +2400,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 @@ -2816,7 +2948,13 @@ c & " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i)) xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) zi=0.5D0*(c(3,i)+c(3,i+1)) - +C Returning the ith atom to box + xi=mod(xi,boxxsize) + if (xi.lt.0) xi=xi+boxxsize + yi=mod(yi,boxysize) + if (yi.lt.0) yi=yi+boxysize + zi=mod(zi,boxzsize) + if (zi.lt.0) zi=zi+boxzsize do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) @@ -2827,28 +2965,73 @@ c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi c zj=c(3,nres+j)-zi C Uncomment following three lines for Ca-p interactions - xj=c(1,j)-xi - yj=c(2,j)-yi - zj=c(3,j)-zi + xj=c(1,j) + yj=c(2,j) + zj=c(3,j) +C returning the jth atom to box + xj=mod(xj,boxxsize) + if (xj.lt.0) xj=xj+boxxsize + yj=mod(yj,boxysize) + if (yj.lt.0) yj=yj+boxysize + zj=mod(zj,boxzsize) + if (zj.lt.0) zj=zj+boxzsize + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 +C Finding the closest jth atom + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif rrij=1.0D0/(xj*xj+yj*yj+zj*zj) +C sss is scaling function for smoothing the cutoff gradient otherwise +C the gradient would not be continuouse + sss=sscale(1.0d0/(dsqrt(rrij))) + if (sss.le.0.0d0) cycle + sssgrad=sscagrad(1.0d0/(dsqrt(rrij))) fac=rrij**expon2 e1=fac*fac*aad(itypj,iteli) e2=fac*bad(itypj,iteli) if (iabs(j-i) .le. 2) then e1=scal14*e1 e2=scal14*e2 - evdw2_14=evdw2_14+e1+e2 + evdw2_14=evdw2_14+(e1+e2)*sss endif evdwij=e1+e2 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 + evdw2=evdw2+evdwij*sss if (calc_grad) then C C Calculate contributions to the gradient in the virtual-bond and SC vectors. C - fac=-(evdwij+e1)*rrij + fac=-(evdwij+e1)*rrij*sss + fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/expon ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac @@ -3083,24 +3266,29 @@ c estr1=0.0d0 c write (iout,*) "distchainmax",distchainmax do i=nnt+1,nct - if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then - estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax) - do j=1,3 - gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) - & *dc(j,i-1)/vbld(i) - enddo - if (energy_dec) write(iout,*) - & "estr1",i,vbld(i),distchainmax, - & gnmr1(vbld(i),-1.0d0,distchainmax) - else + if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle +C estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax) +C do j=1,3 +C gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) +C & *dc(j,i-1)/vbld(i) +C enddo +C if (energy_dec) write(iout,*) +C & "estr1",i,vbld(i),distchainmax, +C & gnmr1(vbld(i),-1.0d0,distchainmax) +C else + if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then + diff = vbld(i)-vbldpDUM + else diff = vbld(i)-vbldp0 c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff + endif estr=estr+diff*diff do j=1,3 gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i) enddo - endif - +C endif +C write (iout,'(a7,i5,4f7.3)') +C & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff enddo estr=0.5d0*AKP*estr+estr1 c @@ -3112,8 +3300,8 @@ c nbi=nbondterm(iti) if (nbi.eq.1) then diff=vbld(i+nres)-vbldsc0(1,iti) -c write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff, -c & AKSC(1,iti),AKSC(1,iti)*diff*diff +C write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff, +C & AKSC(1,iti),AKSC(1,iti)*diff*diff estr=estr+0.5d0*AKSC(1,iti)*diff*diff do j=1,3 gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres) @@ -3142,8 +3330,8 @@ c & AKSC(1,iti),AKSC(1,iti)*diff*diff usum=usum+uprod1 usumsqder=usumsqder+ud(j)*uprod2 enddo -c write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti), -c & AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi) + write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti), + & AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi) estr=estr+uprod/usum do j=1,3 gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres) @@ -3184,7 +3372,10 @@ c write (iout,*) "nres",nres c write (*,'(a,i2)') 'EBEND ICG=',icg c write (iout,*) ithet_start,ithet_end do i=ithet_start,ithet_end - if (itype(i-1).eq.ntyp1) cycle +C if (itype(i-1).eq.ntyp1) cycle + if (i.le.2) cycle + if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 + & .or.itype(i).eq.ntyp1) cycle C Zero the energy function and its derivative at 0 or pi. call splinthet(theta(i),0.5d0*delta,ss,ssd) it=itype(i-1) @@ -3200,8 +3391,12 @@ C Zero the energy function and its derivative at 0 or pi. ichir21=isign(1,itype(i)) ichir22=isign(1,itype(i)) endif + if (i.eq.3) then + y(1)=0.0D0 + y(2)=0.0D0 + else - 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) icrc=0 @@ -3216,7 +3411,8 @@ C Zero the energy function and its derivative at 0 or pi. y(1)=0.0D0 y(2)=0.0D0 endif - if (i.lt.nres .and. itype(i).ne.ntyp1) then + endif + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) icrc=0 @@ -3283,6 +3479,8 @@ C Derivatives of the "mean" values in gamma1 and gamma2. & E_theta,E_tc) endif etheta=etheta+ethetai +c write (iout,'(a6,i5,0pf7.3,f7.3,i5)') +c & 'ebend',i,ethetai,theta(i),itype(i) c write (iout,'(2i3,3f8.3,f10.5)') i,it,rad2deg*theta(i), c & rad2deg*phii,rad2deg*phii1,ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1 @@ -3431,7 +3629,11 @@ 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 (i.eq.2) cycle +C if (itype(i-1).eq.ntyp1) cycle + if (i.le.2) 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 @@ -3443,7 +3645,15 @@ 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.eq.3) then + phii=0.0d0 + ityp1=nthetyp+1 + do k=1,nsingle + cosph1(k)=0.0d0 + sinph1(k)=0.0d0 + enddo + else + if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -3463,7 +3673,8 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) sinph1(k)=0.0d0 enddo endif - if (i.lt.nres .and. itype(i).ne.ntyp1) then + endif + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -3623,14 +3834,14 @@ C ALPHA and OMEGA. common /sccalc/ time11,time12,time112,theti,it,nlobit delta=0.02d0*pi escloc=0.0D0 -c write (iout,'(a)') 'ESC' +C write (iout,*) 'ESC' do i=loc_start,loc_end it=itype(i) if (it.eq.ntyp1) cycle if (it.eq.10) goto 1 nlobit=nlob(iabs(it)) c print *,'i=',i,' it=',it,' nlobit=',nlobit -c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad +C write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad theti=theta(i+1)-pipol x(1)=dtan(theti) x(2)=alph(i) @@ -3666,8 +3877,8 @@ c write (iout,*) "i",i," x",x(1),x(2),x(3) dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k) enddo dersc(2)=dersc(2)+ssd*(escloci-esclocbi) -c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci, -c & esclocbi,ss,ssd + write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci, + & esclocbi,ss,ssd escloci=ss*escloci+(1.0d0-ss)*esclocbi c escloci=esclocbi c write (iout,*) escloci @@ -3701,15 +3912,17 @@ c write (iout,*) escloci enddo dersc(2)=dersc(2)+ssd*(escloci-esclocbi) c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci, -c & esclocbi,ss,ssd +c & esclocbi,ss,ssd escloci=ss*escloci+(1.0d0-ss)*esclocbi -c write (iout,*) escloci +C write (iout,*) 'i=',i, escloci else call enesc(x,escloci,dersc,ddummy,.false.) endif escloc=escloc+escloci -c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc +C write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc + write (iout,'(a6,i5,0pf7.3)') + & 'escloc',i,escloci gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & wscloc*dersc(1) @@ -3946,7 +4159,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) @@ -3986,7 +4199,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 @@ -4029,6 +4242,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 @@ -4158,9 +4373,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)) @@ -4397,8 +4612,11 @@ C Set lprn=.true. for debugging c lprn=.true. etors=0.0D0 do i=iphi_start,iphi_end - if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 - & .or. itype(i).eq.ntyp1) cycle + if (i.le.2) cycle + if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 + & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle +C if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 +C & .or. itype(i).eq.ntyp1) cycle if (itel(i-2).eq.0 .or. itel(i-1).eq.0) goto 1215 if (iabs(itype(i)).eq.20) then iblock=2 @@ -4498,8 +4716,12 @@ C Set lprn=.true. for debugging c lprn=.true. etors_d=0.0D0 do i=iphi_start,iphi_end-1 - if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 - & .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle + if (i.le.3) cycle +C if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 +C & .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle + if ((itype(i-2).eq.ntyp1).or.itype(i-3).eq.ntyp1.or. + & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1).or. + & (itype(i+1).eq.ntyp1)) cycle if (itel(i-2).eq.0 .or. itel(i-1).eq.0 .or. itel(i).eq.0) & goto 1215 itori=itortyp(itype(i-2)) @@ -4731,9 +4953,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 @@ -4756,9 +4978,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 @@ -6497,7 +6720,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 @@ -7507,4 +7730,34 @@ C----------------------------------------------------------------------------- scalar=sc return end +C----------------------------------------------------------------------- + double precision function sscale(r) + double precision r,gamm + include "COMMON.SPLITELE" + if(r.lt.r_cut-rlamb) then + sscale=1.0d0 + else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then + gamm=(r-(r_cut-rlamb))/rlamb + sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0) + else + sscale=0d0 + endif + return + end +C----------------------------------------------------------------------- +C----------------------------------------------------------------------- + double precision function sscagrad(r) + double precision r,gamm + include "COMMON.SPLITELE" + if(r.lt.r_cut-rlamb) then + sscagrad=0.0d0 + else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then + gamm=(r-(r_cut-rlamb))/rlamb + sscagrad=gamm*(6*gamm-6.0d0)/rlamb + else + sscagrad=0.0d0 + endif + return + end +C-----------------------------------------------------------------------