From 0d0c28ca14abcf99ccb778cdba15a9646736571d Mon Sep 17 00:00:00 2001 From: Adam Kazimierz Sieradzan Date: Thu, 15 May 2014 03:53:58 -0400 Subject: [PATCH] zmiany w WHAM PBC przed sprawdzeniem differow --- .gitignore | 2 + source/wham/src-M/COMMON.CHAIN | 4 + source/wham/src-M/COMMON.SPLITELE | 2 + source/wham/src-M/energy_p_new.F | 295 +++++++++++++++++++++++++++++++------ source/wham/src-M/parmread.F | 4 +- source/wham/src-M/readpdb.f | 59 ++++++-- source/wham/src-M/readrtns.F | 7 + 7 files changed, 318 insertions(+), 55 deletions(-) create mode 100644 source/wham/src-M/COMMON.SPLITELE diff --git a/.gitignore b/.gitignore index 2e4973f..d0edd88 100644 --- a/.gitignore +++ b/.gitignore @@ -31,3 +31,5 @@ period_build build_devel build_theta build_new +build_prere/ +period_build2 diff --git a/source/wham/src-M/COMMON.CHAIN b/source/wham/src-M/COMMON.CHAIN index 8dcdd98..b147f08 100644 --- a/source/wham/src-M/COMMON.CHAIN +++ b/source/wham/src-M/COMMON.CHAIN @@ -12,3 +12,7 @@ &nsup,nstart_sup,anatemp, &nend_sup,chain_length,tabperm(maxperm,maxsym),nperm, & nstart_seq,ishift_pdb + double precision boxxsize,boxysize,boxzsize,enecut,sscut,sss, + &sssgrad + common /box/ boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad + diff --git a/source/wham/src-M/COMMON.SPLITELE b/source/wham/src-M/COMMON.SPLITELE new file mode 100644 index 0000000..a2f0447 --- /dev/null +++ b/source/wham/src-M/COMMON.SPLITELE @@ -0,0 +1,2 @@ + double precision r_cut,rlamb + common /splitele/ r_cut,rlamb diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index 8f932cc..a4f5fb4 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -797,6 +797,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 +838,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 +913,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 +944,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 @@ -1852,7 +1905,10 @@ cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e 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 (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 + & .or. itype(i+2).eq.ntyp1 + & .or. itype(i-1).eq.ntyp1 + &) cycle if (itel(i).eq.0) goto 1215 dxi=dc(1,i) dyi=dc(2,i) @@ -1863,10 +1919,21 @@ 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) do j=ielstart(i),ielend(i) - if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle + if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 + & .or.itype(j+2).eq.ntyp1 + & .or.itype(j-1).eq.ntyp1 + &) cycle +C +C) cycle if (itel(j).eq.0) goto 1216 ind=ind+1 iteli=itel(i) @@ -1888,10 +1955,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,7 +2022,7 @@ 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 @@ -1929,7 +2036,7 @@ 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 +2062,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 +2088,7 @@ C enddo enddo #else - facvdw=ev1+evdwij + facvdw=(ev1+evdwij)*sss facel=el1+eesij fac1=fac fac=-3*rrmij*(facvdw+facvdw+facel) @@ -2818,7 +2934,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) @@ -2829,28 +2951,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 @@ -3085,16 +3252,19 @@ 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 estr=estr+diff*diff @@ -3186,7 +3356,9 @@ 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 ((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) @@ -3203,7 +3375,7 @@ C Zero the energy function and its derivative at 0 or pi. ichir22=isign(1,itype(i)) endif - 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 @@ -3218,7 +3390,7 @@ 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 + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) icrc=0 @@ -3433,7 +3605,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 @@ -3445,7 +3619,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 @@ -3465,7 +3639,7 @@ 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 + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -4401,8 +4575,10 @@ 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 (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 @@ -4502,8 +4678,11 @@ 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 +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)) @@ -7512,4 +7691,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----------------------------------------------------------------------- diff --git a/source/wham/src-M/parmread.F b/source/wham/src-M/parmread.F index c9ac2a3..ea6cb14 100644 --- a/source/wham/src-M/parmread.F +++ b/source/wham/src-M/parmread.F @@ -150,7 +150,7 @@ c Read the virtual-bond parameters, masses, and moments of inertia c and Stokes' radii of the peptide group and side chains c #ifdef CRYST_BOND - read (ibond,*) vbldp0,akp + read (ibond,*) vbldp0,vbldpdum,akp do i=1,ntyp nbondterm(i)=1 read (ibond,*) vbldsc0(1,i),aksc(1,i) @@ -162,7 +162,7 @@ c endif enddo #else - read (ibond,*) ijunk,vbldp0,akp,rjunk + read (ibond,*) ijunk,vbldp0,vbldpdum,akp,rjunk do i=1,ntyp read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i), & j=1,nbondterm(i)) diff --git a/source/wham/src-M/readpdb.f b/source/wham/src-M/readpdb.f index a65709e..daddd52 100644 --- a/source/wham/src-M/readpdb.f +++ b/source/wham/src-M/readpdb.f @@ -26,7 +26,9 @@ C geometry. goto 10 else if (card(:3).eq.'TER') then C End current chain - ires_old=ires+1 +c ires_old=ires+1 + ires_old=ires+2 + itype(ires_old-1)=ntyp1 itype(ires_old)=ntyp1 ibeg=2 c write (iout,*) "Chain ended",ires,ishift,ires_old @@ -85,14 +87,51 @@ C system nres=ires do i=2,nres-1 c write (iout,*) i,itype(i) + if (itype(i).eq.ntyp1) then -c write (iout,*) "dummy",i,itype(i) - do j=1,3 - c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2 -c c(j,i)=(c(j,i-1)+c(j,i+1))/2 - dc(j,i)=c(j,i) - enddo - endif + if (itype(i+1).eq.ntyp1) then +C 16/01/2014 by Adasko: Adding to dummy atoms in the chain +C first is connected prevous chain (itype(i+1).eq.ntyp1)=true +C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false +C if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the last dummy residue +C call refsys(i-3,i-2,i-1,e1,e2,e3,fail) +C if (fail) then +C e2(1)=0.0d0 +C e2(2)=1.0d0 +C e2(3)=0.0d0 +C endif !fail +C do j=1,3 +C c(j,i)=c(j,i-1)-1.9d0*e2(j) +C enddo +C else !unres_pdb + do j=1,3 + dcj=(c(j,i-2)-c(j,i-3))/2.0 + c(j,i)=c(j,i-1)+dcj + c(j,nres+i)=c(j,i) + enddo +C endif !unres_pdb + else !itype(i+1).eq.ntyp1 +C if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the first dummy residue +C call refsys(i+1,i+2,i+3,e1,e2,e3,fail) +C if (fail) then +C e2(1)=0.0d0 +C e2(2)=1.0d0 +C e2(3)=0.0d0 +C endif +C do j=1,3 +C c(j,i)=c(j,i+1)-1.9d0*e2(j) +C enddo +C else !unres_pdb + do j=1,3 + dcj=(c(j,i+3)-c(j,i+2))/2.0 + c(j,i)=c(j,i+1)-dcj + c(j,nres+i)=c(j,i) + enddo +C endif !unres_pdb + endif !itype(i+1).eq.ntyp1 + endif !itype.eq.ntyp1 enddo C Calculate the CM of the last side chain. call sccenter(ires,iii,sccor) @@ -102,7 +141,7 @@ C Calculate the CM of the last side chain. nres=nres+1 itype(nres)=ntyp1 do j=1,3 - dcj=c(j,nres-2)-c(j,nres-3) + dcj=(c(j,nres-2)-c(j,nres-3))/2.0 c(j,nres)=c(j,nres-1)+dcj c(j,2*nres)=c(j,nres) enddo @@ -120,7 +159,7 @@ C Calculate the CM of the last side chain. nsup=nsup-1 nstart_sup=2 do j=1,3 - dcj=c(j,4)-c(j,3) + dcj=(c(j,4)-c(j,3))/2.0 c(j,1)=c(j,2)-dcj c(j,nres+1)=c(j,1) enddo diff --git a/source/wham/src-M/readrtns.F b/source/wham/src-M/readrtns.F index 36c13b1..9023170 100644 --- a/source/wham/src-M/readrtns.F +++ b/source/wham/src-M/readrtns.F @@ -17,6 +17,7 @@ include "COMMON.FREE" include "COMMON.CONTROL" include "COMMON.ENERGIES" + include "COMMON.SPLITELE" character*800 controlcard integer i,j,k,ii,n_ene_found integer ind,itype1,itype2,itypf,itypsc,itypp @@ -74,6 +75,12 @@ call readi(controlcard,"RESCALE",rescale_mode,1) check_conf=index(controlcard,"NO_CHECK_CONF").eq.0 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0) + call reada(controlcard,'BOXX',boxxsize,100.0d0) + call reada(controlcard,'BOXY',boxysize,100.0d0) + call reada(controlcard,'BOXZ',boxzsize,100.0d0) +c Cutoff range for interactions + call reada(controlcard,"R_CUT",r_cut,15.0d0) + call reada(controlcard,"LAMBDA",rlamb,0.3d0) call readi(controlcard,'SYM',symetr,1) write (iout,*) "DISTCHAINMAX",distchainmax write (iout,*) "delta",delta -- 1.7.9.5