From: Adam Sieradzan Date: Sat, 8 Feb 2014 17:27:39 +0000 (+0100) Subject: zmiany w galezi multichain X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=0bb81c1c3180a2079d70af7dd534295e1e0b1e4c zmiany w galezi multichain --- diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN index bf5f688..777cc43 100644 --- a/source/unres/src_MD-M/COMMON.CHAIN +++ b/source/unres/src_MD-M/COMMON.CHAIN @@ -14,5 +14,5 @@ & nsup,nstart_sup,nstart_seq, & chain_length,iprzes,tabperm(maxperm,maxsym),nperm common /from_zscore/ nz_start,nz_end,iz_sc - double precision boxxsize,boxysize,boxzsize - common /box/ boxxsize,boxysize,boxzsize + double precision boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad + common /box/ boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad diff --git a/source/unres/src_MD-M/COMMON.DERIV b/source/unres/src_MD-M/COMMON.DERIV index 524d72a..a5e35e2 100644 --- a/source/unres/src_MD-M/COMMON.DERIV +++ b/source/unres/src_MD-M/COMMON.DERIV @@ -9,7 +9,7 @@ & gvdwc(3,maxres),gelc(3,maxres),gelc_long(3,maxres), & gvdwpp(3,maxres),gvdwc_scpp(3,maxres), & gradx_scp(3,maxres),gvdwc_scp(3,maxres),ghpbx(3,maxres), - & ghpbc(3,maxres),gloc(maxvar,2),gradcorr(3,maxres), + & ghpbc(3,maxres),gloc(0:maxvar,2),gradcorr(3,maxres), & gradcorr_long(3,maxres),gradcorr5_long(3,maxres), & gradcorr6_long(3,maxres),gcorr6_turn_long(3,maxres), & gradxorr(3,maxres),gradcorr5(3,maxres),gradcorr6(3,maxres), diff --git a/source/unres/src_MD-M/chainbuild.F b/source/unres/src_MD-M/chainbuild.F index 45a1a53..c4970d3 100644 --- a/source/unres/src_MD-M/chainbuild.F +++ b/source/unres/src_MD-M/chainbuild.F @@ -12,9 +12,121 @@ C include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.INTERACT' - logical lprn + double precision e1(3),e2(3),e3(3) + logical lprn,perbox,fail C Set lprn=.true. for debugging lprn = .false. + perbox=.false. + fail=.false. + if (perbox) then + cost=dcos(theta(3)) + sint=dsin(theta(3)) + print *,'before refsys' + call refsys(2,3,4,e1,e2,e3,fail) + print *,'after refsys' + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + dc(1,0)=c(1,1) + dc(2,0)=c(2,1) + dc(3,0)=c(3,1) + print *,'dc',dc(1,0),dc(2,0),dc(3,0) + dc(1,1)=c(1,2)-c(1,1) + dc(2,1)=c(2,2)-c(2,1) + dc(3,1)=c(3,2)-c(3,1) + dc(1,2)=c(1,3)-c(1,2) + dc(2,2)=c(2,3)-c(2,2) + dc(3,2)=c(3,3)-c(3,2) + t(1,1,1)=e1(1) + t(1,2,1)=e1(2) + t(1,3,1)=e1(3) + t(2,1,1)=e2(1) + t(2,2,1)=e2(2) + t(2,3,1)=e2(3) + t(3,1,1)=e3(1) + t(3,2,1)=e3(2) + t(3,3,1)=e3(3) + veclen=0.0d0 + do i=1,3 + veclen=veclen+(c(i,2)-c(i,1))**2 + enddo + veclen=sqrt(veclen) + r(1,1,1)= 1.0D0 + r(1,2,1)= 0.0D0 + r(1,3,1)= 0.0D0 + r(2,1,1)= 0.0D0 + r(2,2,1)= 1.0D0 + r(2,3,1)= 0.0D0 + r(3,1,1)= 0.0D0 + r(3,2,1)= 0.0D0 + r(3,3,1)= 1.0D0 + do i=1,3 + do j=1,3 + rt(i,j,1)=t(i,j,1) + enddo + enddo + do i=1,3 + do j=1,3 + prod(i,j,1)=0.0D0 + prod(i,j,2)=t(i,j,1) + enddo + prod(i,i,1)=1.0D0 + enddo + call locate_side_chain(2) + do i=4,nres +#ifdef OSF + theti=theta(i) + if (theti.ne.theti) theti=100.0 + phii=phi(i) + if (phii.ne.phii) phii=180.0 +#else + theti=theta(i) + phii=phi(i) +#endif + cost=dcos(theti) + sint=dsin(theti) + cosphi=dcos(phii) + sinphi=dsin(phii) +* Define the matrices of the rotation about the virtual-bond valence angles +* theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this +* program), R(i,j,k), and, the cumulative matrices of rotation RT + t(1,1,i-2)=-cost + t(1,2,i-2)=-sint + t(1,3,i-2)= 0.0D0 + t(2,1,i-2)=-sint + t(2,2,i-2)= cost + t(2,3,i-2)= 0.0D0 + t(3,1,i-2)= 0.0D0 + t(3,2,i-2)= 0.0D0 + t(3,3,i-2)= 1.0D0 + r(1,1,i-2)= 1.0D0 + r(1,2,i-2)= 0.0D0 + r(1,3,i-2)= 0.0D0 + r(2,1,i-2)= 0.0D0 + r(2,2,i-2)=-cosphi + r(2,3,i-2)= sinphi + r(3,1,i-2)= 0.0D0 + r(3,2,i-2)= sinphi + r(3,3,i-2)= cosphi + rt(1,1,i-2)=-cost + rt(1,2,i-2)=-sint + rt(1,3,i-2)=0.0D0 + rt(2,1,i-2)=sint*cosphi + rt(2,2,i-2)=-cost*cosphi + rt(2,3,i-2)=sinphi + rt(3,1,i-2)=-sint*sinphi + rt(3,2,i-2)=cost*sinphi + rt(3,3,i-2)=cosphi + call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1)) + do j=1,3 + dc_norm(j,i-1)=prod(j,1,i-1) + dc(j,i-1)=vbld(i)*prod(j,1,i-1) + enddo + call locate_side_chain(i-1) + enddo + else C C Define the origin and orientation of the coordinate system and locate the C first three CA's and SC(2). @@ -59,7 +171,7 @@ C 1212 format (a3,'(',i3,')',2(f10.5,2f10.2)) endif - + endif return end c------------------------------------------------------------------------- @@ -126,8 +238,8 @@ C dc_norm(3,1)=0.0D0 do j=1,3 dc_norm(j,2)=prod(j,1,2) - dc(j,2)=vbld(3)*prod(j,1,2) - c(j,3)=c(j,2)+dc(j,2) + dc(j,2)=vbld(3)*prod(j,1,2) + c(j,3)=c(j,2)+dc(j,2) enddo call locate_side_chain(2) return diff --git a/source/unres/src_MD-M/checkder_p.F b/source/unres/src_MD-M/checkder_p.F index 0539e48..32d2366 100644 --- a/source/unres/src_MD-M/checkder_p.F +++ b/source/unres/src_MD-M/checkder_p.F @@ -272,8 +272,8 @@ C Check the gradient of the energy in Cartesian coordinates. integer uiparm(1) double precision urparm(1) external fdum - r_cut=2.0d0 - rlambd=0.3d0 +c r_cut=2.0d0 +c rlambd=0.3d0 icg=1 nf=0 nfl=0 diff --git a/source/unres/src_MD-M/contact.f b/source/unres/src_MD-M/contact.f index 24b11d6..cc4e0b7 100644 --- a/source/unres/src_MD-M/contact.f +++ b/source/unres/src_MD-M/contact.f @@ -175,7 +175,7 @@ c do i=1,nharp c write (iout,*)'i',i,' iharp',(iharp(k,i),k=1,4) c enddo if (lprint) then - write (iout,*) "Hairpins:" + write (iout,*) "Hairpins:",nharp do i=1,nharp i1=iharp(1,i) j1=iharp(2,i) diff --git a/source/unres/src_MD-M/elecont.f b/source/unres/src_MD-M/elecont.f index a962630..4faf3e5 100644 --- a/source/unres/src_MD-M/elecont.f +++ b/source/unres/src_MD-M/elecont.f @@ -214,7 +214,7 @@ c-------------------------------------------- double precision p1,p2 external freeres - if(.not.dccart) call chainbuild +cc???? if(.not.dccart) call chainbuild cd call write_pdb(99,'sec structure',0d0) ncont=0 nbfrag=0 diff --git a/source/unres/src_MD-M/energy_p_new-sep.F b/source/unres/src_MD-M/energy_p_new-sep.F index 0b8f27b..85f3e1f 100644 --- a/source/unres/src_MD-M/energy_p_new-sep.F +++ b/source/unres/src_MD-M/energy_p_new-sep.F @@ -2,6 +2,7 @@ C----------------------------------------------------------------------- double precision function sscale(r) double precision r,gamm include "COMMON.SPLITELE" + include "COMMON.CHAIN" if(r.lt.r_cut-rlamb) then sscale=1.0d0 else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then @@ -13,6 +14,23 @@ C----------------------------------------------------------------------- return end C----------------------------------------------------------------------- +C----------------------------------------------------------------------- + double precision function sscagrad(r) + double precision r,gamm + include "COMMON.SPLITELE" + include "COMMON.CHAIN" + 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----------------------------------------------------------------------- + subroutine elj_long(evdw) C C This subroutine calculates the interaction energy of nonbonded side chains diff --git a/source/unres/src_MD-M/energy_p_new-sep_barrier.F b/source/unres/src_MD-M/energy_p_new-sep_barrier.F index 6592ace..574e19b 100644 --- a/source/unres/src_MD-M/energy_p_new-sep_barrier.F +++ b/source/unres/src_MD-M/energy_p_new-sep_barrier.F @@ -13,6 +13,21 @@ C----------------------------------------------------------------------- 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----------------------------------------------------------------------- subroutine elj_long(evdw) C C This subroutine calculates the interaction energy of nonbonded side chains 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 635a41e..8f89ae9 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -24,6 +24,7 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.MD' include 'COMMON.CONTROL' include 'COMMON.TIME1' + include 'COMMON.SPLITELE' #ifdef MPI c print*,"ETOTAL Processor",fg_rank," absolute rank",myrank, c & " nfgtasks",nfgtasks @@ -1412,6 +1413,7 @@ C include 'COMMON.IOUNITS' include 'COMMON.CALC' include 'COMMON.CONTROL' + include 'COMMON.SPLITELE' logical lprn integer xshift,yshift,zshift evdw=0.0D0 @@ -1451,8 +1453,8 @@ C Condition for being inside the proper box go to 135 endif 136 continue - if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxxsize - if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxxsize + if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize + if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize C Condition for being inside the proper box if ((zi.gt.((zshift+0.5d0)*boxzsize)).or. & (zi.lt.((zshift-0.5d0)*boxzsize))) then @@ -1520,8 +1522,8 @@ C Condition for being inside the proper box go to 138 endif 139 continue - if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize - if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize + if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize + if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize C Condition for being inside the proper box if ((zj.gt.((0.5d0)*boxzsize)).or. & (zj.lt.((-0.5d0)*boxzsize))) then @@ -1539,6 +1541,12 @@ c write (iout,*) "j",j," dc_norm", c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j) 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)) + +c write (iout,'(a7,4f8.3)') +c & "ssscale",sss,((1.0d0/rij)/sigma(itypi,itypj)),r_cut,rlamb + if (sss.gt.0.0d0) then C Calculate angle-dependent terms of energy and contributions to their C derivatives. call sc_angular @@ -1567,7 +1575,7 @@ c--------------------------------------------------------------- c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt, c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 evdwij=evdwij*eps2rt*eps3rt - evdw=evdw+evdwij + evdw=evdw+evdwij*sss if (lprn) then sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) epsi=bb(itypi,itypj)**2/aa(itypi,itypj) @@ -1587,6 +1595,9 @@ C Calculate gradient components. fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder fac=rij*fac +c print '(2i4,6f8.4)',i,j,sss,sssgrad* +c & evdwij,fac,sigma(itypi,itypj),expon + fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij c fac=0.0d0 C Calculate the radial part of the gradient gg(1)=xj*fac @@ -1594,6 +1605,7 @@ C Calculate the radial part of the gradient gg(3)=zj*fac C Calculate angular part of the gradient. call sc_grad + endif enddo ! j enddo ! iint enddo ! i @@ -1807,6 +1819,7 @@ C---------------------------------------------------------------------------- include 'COMMON.CALC' include 'COMMON.IOUNITS' double precision dcosom1(3),dcosom2(3) +cc print *,'sss=',sss eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 @@ -1825,16 +1838,16 @@ c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k)) enddo do k=1,3 - gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k) + gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*sss enddo c write (iout,*) "gg",(gg(k),k=1,3) do k=1,3 gvdwx(k,i)=gvdwx(k,i)-gg(k) & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) - & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv + & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*sss gvdwx(k,j)=gvdwx(k,j)+gg(k) & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) - & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv + & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*sss c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) @@ -2770,6 +2783,7 @@ C include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.TIME1' + include 'COMMON.SPLITELE' dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), @@ -2883,8 +2897,8 @@ C Condition for being inside the proper box go to 185 endif 186 continue - if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxxsize - if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxxsize + if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize + if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize C Condition for being inside the proper box if ((zmedi.gt.((0.5d0)*boxzsize)).or. & (zmedi.lt.((-0.5d0)*boxzsize))) then @@ -2926,8 +2940,8 @@ C Condition for being inside the proper box go to 195 endif 196 continue - if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxxsize - if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxxsize + if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize + if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize C Condition for being inside the proper box if ((zmedi.gt.((0.5d0)*boxzsize)).or. & (zmedi.lt.((-0.5d0)*boxzsize))) then @@ -2976,8 +2990,8 @@ C Condition for being inside the proper box go to 165 endif 166 continue - if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxxsize - if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxxsize + if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxzsize + if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize C Condition for being inside the proper box if ((zmedi.gt.((zshift+0.5d0)*boxzsize)).or. & (zmedi.lt.((zshift-0.5d0)*boxzsize))) then @@ -3027,6 +3041,7 @@ C------------------------------------------------------------------------------- include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.TIME1' + include 'COMMON.SPLITELE' dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), @@ -3085,8 +3100,8 @@ C Condition for being inside the proper box go to 175 endif 176 continue - if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize - if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize + if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize + if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize C Condition for being inside the proper box if ((zj.gt.((0.5d0)*boxzsize)).or. & (zj.lt.((-0.5d0)*boxzsize))) then @@ -3097,6 +3112,10 @@ C endif !endPBC condintion yj=yj-ymedi zj=zj-zmedi rij=xj*xj+yj*yj+zj*zj + + sss=sscale(sqrt(rij)) + sssgrad=sscagrad(sqrt(rij)) +c if (sss.gt.0.0d0) then rrmij=1.0D0/rij rij=dsqrt(rij) rmij=1.0D0/rij @@ -3112,14 +3131,15 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions ev2=bbb*r6ij fac3=ael6i*r6ij fac4=ael3i*r3ij - evdwij=ev1+ev2 + evdwij=(ev1+ev2) el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)) el2=fac4*fac - eesij=el1+el2 +C MARYSIA + eesij=(el1+el2) 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 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, @@ -3136,7 +3156,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 @@ -3188,10 +3208,11 @@ cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l) cgrad enddo cgrad enddo #else - facvdw=ev1+evdwij - facel=el1+eesij +C MARYSIA + facvdw=(ev1+evdwij)*sss + facel=(el1+eesij) fac1=fac - fac=-3*rrmij*(facvdw+facvdw+facel) + fac=-3*rrmij*(facvdw+facvdw+facel)+sssgrad*rmij*evdwij erij(1)=xj*rmij erij(2)=yj*rmij erij(3)=zj*rmij @@ -3220,9 +3241,9 @@ cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo c 9/28/08 AL Gradient compotents will be summed only at the end - ggg(1)=facvdw*xj - ggg(2)=facvdw*yj - ggg(3)=facvdw*zj + ggg(1)=facvdw*xj*sss + ggg(2)=facvdw*yj*sss + ggg(3)=facvdw*zj*sss do k=1,3 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) @@ -3261,14 +3282,16 @@ cgrad enddo cgrad enddo do k=1,3 gelc(k,i)=gelc(k,i) - & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) - & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) + & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) + & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) gelc(k,j)=gelc(k,j) - & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) - & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) + & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) + & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo +C MARYSIA +c endif !sscale IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN @@ -3460,9 +3483,9 @@ cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eelloc',i,j,eel_loc_ij - if (eel_loc_ij.ne.0) - & write (iout,'(a4,2i4,8f9.5)')'chuj', - & i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4) +c if (eel_loc_ij.ne.0) +c & write (iout,'(a4,2i4,8f9.5)')'chuj', +c & i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4) eel_loc=eel_loc+eel_loc_ij C Partial derivatives in virtual-bond dihedral angles gamma @@ -3490,14 +3513,14 @@ cgrad enddo cgrad enddo C Remaining derivatives of eello do l=1,3 - gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+ - & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4) - gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+ - & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4) - gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+ - & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4) - gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+ - & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4) + gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+ + & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)) + gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+ + & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)) + gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+ + & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)) + gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+ + & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)) enddo ENDIF C Change 12/26/95 to calculate four-body contributions to H-bonding energy @@ -4050,8 +4073,8 @@ C Condition for being inside the proper box go to 135 endif 136 continue - if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxxsize - if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxxsize + if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize + if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize C Condition for being inside the proper box if ((zi.gt.((zshift+0.5d0)*boxzsize)).or. & (zi.lt.((zshift-0.5d0)*boxzsize))) then @@ -4087,8 +4110,8 @@ C Condition for being inside the proper box go to 175 endif 176 continue - if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize - if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize + if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize + if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize C Condition for being inside the proper box if ((zj.gt.((0.5d0)*boxzsize)).or. & (zj.lt.((-0.5d0)*boxzsize))) then @@ -4098,6 +4121,7 @@ C Condition for being inside the proper box yj=yj-yi zj=zj-zi rij=xj*xj+yj*yj+zj*zj + r0ij=r0_scp r0ijsq=r0ij*r0ij if (rij.lt.r0ijsq) then @@ -4171,6 +4195,7 @@ C include 'COMMON.FFIELD' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' + include 'COMMON.SPLITELE' dimension ggg(3) evdw2=0.0D0 evdw2_14=0.0d0 @@ -4203,8 +4228,8 @@ C Condition for being inside the proper box go to 135 endif 136 continue - if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxxsize - if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxxsize + if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize + if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize C Condition for being inside the proper box if ((zi.gt.((zshift+0.5d0)*boxzsize)).or. & (zi.lt.((zshift-0.5d0)*boxzsize))) then @@ -4240,8 +4265,8 @@ C Condition for being inside the proper box go to 175 endif 176 continue - if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize - if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize + if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize + if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize C Condition for being inside the proper box if ((zj.gt.((0.5d0)*boxzsize)).or. & (zj.lt.((-0.5d0)*boxzsize))) then @@ -4251,23 +4276,27 @@ C Condition for being inside the proper box yj=yj-yi zj=zj-zi rrij=1.0D0/(xj*xj+yj*yj+zj*zj) + sss=sscale(1.0d0/(dsqrt(rrij))) + sssgrad=sscagrad(1.0d0/(dsqrt(rrij))) + if (sss.gt.0.0d0) then 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 - evdw2=evdw2+evdwij + evdw2=evdw2+evdwij*sss 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 - 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 @@ -4302,7 +4331,8 @@ cgrad enddo gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k) gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) enddo - enddo + endif !endif for sscale cutoff + enddo ! j enddo ! iint enddo ! i @@ -4622,7 +4652,9 @@ c time12=1.0d0 etheta=0.0D0 c write (*,'(a,i2)') 'EBEND ICG=',icg do i=ithet_start,ithet_end + print *,i,itype(i-1),itype(i),itype(i-2) if (itype(i-1).eq.ntyp1) cycle + print *,'wchodze',itype(i-1) 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) @@ -4720,7 +4752,7 @@ C Derivatives of the "mean" values in gamma1 and gamma2. & 'ebend',i,ethetai,theta(i),itype(i) 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) + gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg) enddo C Ufff.... We've done all this!!! return @@ -4739,7 +4771,7 @@ C Calculate the contributions to both Gaussian lobes. C 6/6/97 - Deform the Gaussians using the factor of 1/(1+time) C The "polynomial part" of the "standard deviation" of this part of C the distributioni. - write (iout,*) thetai,thet_pred_mean +ccc write (iout,*) thetai,thet_pred_mean sig=polthet(3,it) do j=2,0,-1 sig=sig*thet_pred_mean+polthet(j,it) @@ -4864,6 +4896,7 @@ C logical lprn /.false./, lprn1 /.false./ etheta=0.0D0 do i=ithet_start,ithet_end +c print *,i,itype(i-1),itype(i),itype(i-2) 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 @@ -5029,7 +5062,7 @@ 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 - gloc(nphi+i-2,icg)=wang*dethetai + gloc(nphi+i-2,icg)=wang*dethetai+ gloc(nphi+i-2,icg) enddo return end @@ -5892,9 +5925,11 @@ c lprn=.true. etors=0.0D0 do i=iphi_start,iphi_end C ANY TWO ARE DUMMY ATOMS in row CYCLE - if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or. - & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)) .or. - & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle +c if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or. +c & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)) .or. +c & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle + if ((itype(i-3).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. + & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1)) cycle C For introducing the NH3+ and COO- group please check the etor_d for reference C and guidance etors_ii=0.0D0 diff --git a/source/unres/src_MD-M/geomout.F b/source/unres/src_MD-M/geomout.F index c6a3944..4673c4e 100644 --- a/source/unres/src_MD-M/geomout.F +++ b/source/unres/src_MD-M/geomout.F @@ -91,7 +91,7 @@ cmodel write (iunit,'(a5,i6)') 'MODEL',1 ires=0 do i=nnt,nct iti=itype(i) - if (iti.eq.ntyp1) then + if ((iti.eq.ntyp1).and.((itype(i+1)).eq.ntyp1)) then ichain=ichain+1 ires=0 write (iunit,'(a)') 'TER' @@ -99,6 +99,7 @@ cmodel write (iunit,'(a5,i6)') 'MODEL',1 ires=ires+1 iatom=iatom+1 ica(i)=iatom + if (iti.ne.ntyp1) then write (iunit,10) iatom,restyp(iti),chainid(ichain), & ires,(c(j,i),j=1,3),vtot(i) if (iti.ne.10) then @@ -106,6 +107,7 @@ cmodel write (iunit,'(a5,i6)') 'MODEL',1 write (iunit,20) iatom,restyp(iti),chainid(ichain), & ires,(c(j,nres+i),j=1,3), & vtot(i+nres) + endif endif endif enddo diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F index a650241..70e35ab 100644 --- a/source/unres/src_MD-M/initialize_p.F +++ b/source/unres/src_MD-M/initialize_p.F @@ -267,8 +267,8 @@ C C Initialize constants used to split the energy into long- and short-range C components C - r_cut=2.0d0 - rlamb=0.3d0 +C r_cut=2.0d0 +C rlamb=0.3d0 #ifndef SPLITELE nprint_ene=nprint_ene-1 #endif diff --git a/source/unres/src_MD-M/intcartderiv.F b/source/unres/src_MD-M/intcartderiv.F index 369a4f0..d715c19 100644 --- a/source/unres/src_MD-M/intcartderiv.F +++ b/source/unres/src_MD-M/intcartderiv.F @@ -692,7 +692,7 @@ c Check omega gradient enddo return end - +c------------------------------------------------------------ subroutine chainbuild_cart implicit real*8 (a-h,o-z) include 'DIMENSIONS' @@ -736,6 +736,7 @@ c call flush(iout) #endif do j=1,3 c(j,1)=dc(j,0) +c c(j,1)=c(j,1) enddo do i=2,nres do j=1,3 diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index bd2165b..22236bc 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -591,6 +591,7 @@ C read (itorp,*,end=113,err=113) ntortyp,nterm_old if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp) + itortyp(ntyp1)=0 do i=1,ntortyp do j=1,ntortyp read (itorp,'(a)') diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index de4b92e..8e7b152 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -8,6 +8,7 @@ include 'COMMON.CONTROL' include 'COMMON.SBRIDGE' include 'COMMON.IOUNITS' + include 'COMMON.SPLITELE' logical file_exist C Read force-field parameters except weights call parmread @@ -79,6 +80,7 @@ C include 'COMMON.FFIELD' include 'COMMON.INTERACT' include 'COMMON.SETUP' + include 'COMMON.SPLITELE' COMMON /MACHSW/ KDIAG,ICORFL,IXDR character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/ character*80 ucase @@ -219,7 +221,9 @@ C Reading the dimensions of box in x,y,z coordinates 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) if (me.eq.king .or. .not.out1file ) & write (iout,*) "DISTCHAINMAX",distchainmax @@ -346,8 +350,8 @@ C ntime_split0=ntime_split call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64) ntime_split0=ntime_split - call reada(controlcard,"R_CUT",r_cut,2.0d0) - call reada(controlcard,"LAMBDA",rlamb,0.3d0) +c call reada(controlcard,"R_CUT",r_cut,2.0d0) +c call reada(controlcard,"LAMBDA",rlamb,0.3d0) rest = index(controlcard,"REST").gt.0 tbf = index(controlcard,"TBF").gt.0 usampl = index(controlcard,"USAMPL").gt.0 diff --git a/source/unres/src_MD-M/refsys.f b/source/unres/src_MD-M/refsys.f index 86b0524..dff4aa7 100644 --- a/source/unres/src_MD-M/refsys.f +++ b/source/unres/src_MD-M/refsys.f @@ -14,12 +14,16 @@ c form a linear fragment. Returns vectors e1, e2, and e3. include 'COMMON.IOUNITS' include 'COMMON.CHAIN' double precision coinc/1.0D-13/,align /1.0D-13/ +c print *,'just initialize' fail=.false. +c print *,fail s1=0.0 s2=0.0 + print *,s1,s2 do 1 i=1,3 zi=c(i,i2)-c(i,i3) ui=c(i,i4)-c(i,i3) + print *,zi,ui s1=s1+zi*zi s2=s2+ui*ui z(i)=zi @@ -38,6 +42,7 @@ c 3 c(i,i1)=0.0D0 do 5 i=1,3 5 c(i,i1)=0.0D0 return + print *,'two if pass' 4 s1=1.0/s1 s2=1.0/s2 v1=z(2)*u(3)-z(3)*u(2) @@ -60,6 +65,7 @@ c 7 c(i,i1)=0.0D0 e2(1)=e1(3)*e3(2)-e1(2)*e3(3) e2(2)=e1(1)*e3(3)-e1(3)*e3(1) e2(3)=e1(2)*e3(1)-e1(1)*e3(2) + print *,'just before leave' 1000 format (/1x,' * * * Error - atoms',i4,' and',i4,' coincide.', 1 'coordinates of atom',i4,' are set to zero.') 1010 format (/1x,' * * * Error - atoms',2(i4,2h, ),i4,' form a linear', diff --git a/source/unres/src_MD-M/unres.F b/source/unres/src_MD-M/unres.F index 0039fcc..83afd9f 100644 --- a/source/unres/src_MD-M/unres.F +++ b/source/unres/src_MD-M/unres.F @@ -189,9 +189,16 @@ c--------------------------------------------------------------------------- double precision energy(0:n_ene) double precision energy_long(0:n_ene),energy_short(0:n_ene) double precision varia(maxvar) - if (indpdb.eq.0) call chainbuild + if (indpdb.eq.0) call chainbuild time00=MPI_Wtime() + print *,'dc',c(1,1) + if (indpdb.ne.0) then + dc(1,0)=c(1,1) + dc(2,0)=c(2,1) + dc(3,0)=c(3,1) + endif call chainbuild_cart + print *,'dc',dc(1,0),dc(2,0),dc(3,0) if (split_ene) then print *,"Processor",myrank," after chainbuild" icall=1 @@ -216,7 +223,9 @@ c--------------------------------------------------------------------------- etot =etota call enerprint(energy(0)) call hairpin(.true.,nharp,iharp) + print *,'after hairpin' call secondary2(.true.) + print *,'after secondary' if (minim) then crc overlap test if (overlapsc) then @@ -248,8 +257,10 @@ crc overlap test evals=nfun/(MPI_WTIME()-time1) print *,'# eval/s',evals print *,'refstr=',refstr - call hairpin(.true.,nharp,iharp) + call hairpin(.false.,nharp,iharp) + print *,'after hairpin' call secondary2(.true.) + print *,'after secondary' call etotal(energy(0)) etot = energy(0) call enerprint(energy(0))