From bf8b3ab324569ca7ae7bb9045ad9358f6826a150 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Mon, 6 Jul 2015 11:17:47 +0200 Subject: [PATCH] intoduction of quartic restrains in multichain, bugfix in single chain --- PARAM/bond_AM1_ext.parm | 1 - PARAM/pot_theta_G631_DIL.parm | 2 +- source/unres/src_MD-M/checkder_p.F | 12 ++++++++++-- source/unres/src_MD-M/energy_p_new_barrier.F | 13 +++++++++---- source/unres/src_MD-M/gnmr1.f | 4 ++-- source/unres/src_MD-M/intcor.f | 4 ++++ source/unres/src_MD-M/parmread.F | 5 +++-- source/unres/src_MD-M/readpdb.F | 2 ++ source/unres/src_MD-M/readrtns_CSA.F | 18 ++++++++++++++---- source/unres/src_MD-M/unres.F | 6 ++++++ source/unres/src_MD/COMMON.TORSION | 14 +++++++++----- source/unres/src_MD/parmread.F | 25 ++++++++++++++----------- source/unres/src_MD/ssMD.F | 2 +- 13 files changed, 75 insertions(+), 33 deletions(-) diff --git a/PARAM/bond_AM1_ext.parm b/PARAM/bond_AM1_ext.parm index 333a43f..0f7bb5e 100644 --- a/PARAM/bond_AM1_ext.parm +++ b/PARAM/bond_AM1_ext.parm @@ -23,4 +23,3 @@ 1 3.799 124.6 0.000 149.0 49.67 7.2 ! Dap(Bz) 1 0.743 353.00 0.000 42.0 14.00 4.7 ! Aib 1 1.210 353.00 0.000 42.0 14.00 5.6 ! Abu - diff --git a/PARAM/pot_theta_G631_DIL.parm b/PARAM/pot_theta_G631_DIL.parm index 3b94e8b..494bff1 100644 --- a/PARAM/pot_theta_G631_DIL.parm +++ b/PARAM/pot_theta_G631_DIL.parm @@ -1,5 +1,5 @@ 2 10 4 4 6 4 -1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 2 +1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 2 3 Gppreg 4.58415E+01 1.05721E+03 diff --git a/source/unres/src_MD-M/checkder_p.F b/source/unres/src_MD-M/checkder_p.F index b017e1c..4ebfd05 100644 --- a/source/unres/src_MD-M/checkder_p.F +++ b/source/unres/src_MD-M/checkder_p.F @@ -277,11 +277,12 @@ C Check the gradient of the energy in Cartesian coordinates. icg=1 nf=0 nfl=0 + print *,"ATU 3" call intout c call intcartderiv c call checkintcartgrad call zerograd - aincr=10.0D-6 + aincr=8.0D-6 write(iout,*) 'Calling CHECK_ECARTINT.' nf=0 icall=0 @@ -507,8 +508,10 @@ c------------------------------------------------------------------------- #else do i=2,nres #endif +C print *,i dnorm1=dist(i-1,i) - dnorm2=dist(i,i+1) + dnorm2=dist(i,i+1) +C print *,i,dnorm1,dnorm2 do j=1,3 c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1 & +(c(j,i+1)-c(j,i))/dnorm2) @@ -529,11 +532,16 @@ c------------------------------------------------------------------------- endif endif omeg(i)=beta(nres+i,i,maxres2,i+1) +C print *,omeg(i) alph(i)=alpha(nres+i,i,maxres2) +C print *,alph(i) theta(i+1)=alpha(i-1,i,i+1) vbld(i)=dist(i-1,i) +C print *,vbld(i) vbld_inv(i)=1.0d0/vbld(i) vbld(nres+i)=dist(nres+i,i) +C print *,vbld(i+nres) + if (itype(i).ne.10) then vbld_inv(nres+i)=1.0d0/vbld(nres+i) else 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 cf9ddcb..745f7d5 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -4087,6 +4087,9 @@ C include 'COMMON.CONTROL' dimension ggg(3) ehpb=0.0D0 + do i=1,3 + ggg(i)=0.0d0 + enddo cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr cd write(iout,*)'link_start=',link_start,' link_end=',link_end if (link_end.eq.0) return @@ -4122,9 +4125,9 @@ cd & ' waga=',waga,' fac=',fac c Restraints from contact prediction dd=dist(ii,jj) if (constr_dist.eq.11) then - ehpb=ehpb+fordepth(i)**4 + ehpb=ehpb+fordepth(i)**4.0d0 & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) - fac=fordepth(i)**4 + fac=fordepth(i)**4.0d0 & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd else if (dhpb1(i).gt.0.0d0) then @@ -4162,8 +4165,10 @@ C Calculate the distance between the two points and its difference from the C target distance. dd=dist(ii,jj) if (constr_dist.eq.11) then - ehpb=ehpb+fordepth(i)**4*rlornmr1(dd,dhpb(i),dhpb1(i)) - fac=fordepth(i)**4*rlornmr1prim(dd,dhpb(i),dhpb1(i))/dd + ehpb=ehpb+fordepth(i)**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd else if (dhpb1(i).gt.0.0d0) then ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) diff --git a/source/unres/src_MD-M/gnmr1.f b/source/unres/src_MD-M/gnmr1.f index f44e404..8bfc43a 100644 --- a/source/unres/src_MD-M/gnmr1.f +++ b/source/unres/src_MD-M/gnmr1.f @@ -61,10 +61,10 @@ c------------------------------------------------------------------------------ double precision wykl /4.0d0/ if (y.lt.ymin) then rlornmr1prim=-(ymin-y)**(wykl-1)*sigma**wykl*wykl/ - & ((ymin-y)**wykl+sigma**wykl) + & ((ymin-y)**wykl+sigma**wykl)**2 else if (y.gt.ymax) then rlornmr1prim=(y-ymax)**(wykl-1)*sigma**wykl*wykl/ - & ((y-ymax)**wykl+sigma**wykl) + & ((y-ymax)**wykl+sigma**wykl)**2 else rlornmr1prim=0.0d0 endif diff --git a/source/unres/src_MD-M/intcor.f b/source/unres/src_MD-M/intcor.f index a3cd5d0..9195f8a 100644 --- a/source/unres/src_MD-M/intcor.f +++ b/source/unres/src_MD-M/intcor.f @@ -17,7 +17,11 @@ c z23=c(3,i3)-c(3,i2) vnorm=dsqrt(x12*x12+y12*y12+z12*z12) wnorm=dsqrt(x23*x23+y23*y23+z23*z23) + if ((vnorm.eq.0.0).or.(wnorm.eq.0.0)) then + scalar=1.0 + else scalar=(x12*x23+y12*y23+z12*z23)/(vnorm*wnorm) + endif alpha=arcos(scalar) return end diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 7626334..3e9a516 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -73,6 +73,7 @@ c #else read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok do i=1,ntyp + print *,i read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i), & j=1,nbondterm(i)),msc(i),isc(i),restok(i) dsc(i) = vbldsc0(1,i) @@ -924,8 +925,8 @@ c b1(2,i)=0.0d0 B1tilde(2,i) =-b(5) B1tilde(1,-i) =-b(3) B1tilde(2,-i) =b(5) -c b1tilde(1,i)=0.0d0 -c b1tilde(2,i)=0.0d0 + b1tilde(1,i)=0.0d0 + b1tilde(2,i)=0.0d0 B2(1,i) = b(2) B2(2,i) = b(4) B2(1,-i) =b(2) diff --git a/source/unres/src_MD-M/readpdb.F b/source/unres/src_MD-M/readpdb.F index f2daadd..71daf6f 100644 --- a/source/unres/src_MD-M/readpdb.F +++ b/source/unres/src_MD-M/readpdb.F @@ -199,6 +199,7 @@ C Calculate internal coordinates. & (c(j,nres+ires),j=1,3) enddo endif +C print *,"before int_from_cart" call int_from_cart(.true.,.false.) call sc_loc_geom(.true.) do i=1,nres @@ -390,6 +391,7 @@ c vbld(nres)=3.8d0 c vbld_inv(nres)=1.0d0/vbld(2) c endif c endif + print *,"A TU2" if (lside) then do i=2,nres-1 do j=1,3 diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index f53da12..00e72a0 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -708,7 +708,11 @@ C 12/1/95 Added weight for the multi-body term WCORR v2ss=v2ss*wstrain/wsc v3ss=v3ss*wstrain/wsc else - ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain + if (wstrain.ne.0.0) then + ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain + else + ss_depth=0.0 + endif endif if(me.eq.king.or..not.out1file) then @@ -931,7 +935,9 @@ czscore call geom_to_var(nvar,coord_exp_zs(1,1)) enddo call contact(.true.,ncont_ref,icont_ref,co) endif -c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup + endif + print *, "A TU" + write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup call flush(iout) if (constr_dist.gt.0) call read_dist_constr write (iout,*) "After read_dist_constr nhpb",nhpb @@ -951,7 +957,7 @@ c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i) enddo endif - endif +C endif if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and. & modecalc.ne.10) then @@ -1134,6 +1140,7 @@ cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar) & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') & 'Processor',myrank,': end reading molecular data.' #endif + print *,"A TU?" return end c-------------------------------------------------------------------------- @@ -2257,7 +2264,8 @@ c------------------------------------------------------------------------------- integer ifrag_(2,100),ipair_(2,100) double precision wfrag_(100),wpair_(100) character*500 controlcard -c write (iout,*) "Calling read_dist_constr" + print *, "WCHODZE" + write (iout,*) "Calling read_dist_constr" c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup c call flush(iout) call card_concat(controlcard) @@ -2351,12 +2359,14 @@ c write (iout,*) "j",j," k",k enddo endif enddo + print *,ndist_ do i=1,ndist_ if (constr_dist.eq.11) then read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1) fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) else +C print *,"in else" read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), & ibecarb(i),forcon(nhpb+1) endif diff --git a/source/unres/src_MD-M/unres.F b/source/unres/src_MD-M/unres.F index 537a965..ee926b3 100644 --- a/source/unres/src_MD-M/unres.F +++ b/source/unres/src_MD-M/unres.F @@ -676,6 +676,7 @@ c--------------------------------------------------------------------------- include 'COMMON.SBRIDGE' common /srutu/ icall double precision energy(0:max_ene) + print *,"A TU?" c do i=2,nres c vbld(i)=vbld(i)+ran_number(-0.1d0,0.1d0) c if (itype(i).ne.10) @@ -701,10 +702,14 @@ c enddo totT=1.d0 eq_time=0.0d0 call read_fragments + print *, "AFTER read fragments" call chainbuild_cart + print *,"chainbuild_cart" call cartprint + print *,"After cartprint" call intout icall=1 + print *,"before ETOT" call etotal(energy(0)) etot = energy(0) call enerprint(energy(0)) @@ -712,6 +717,7 @@ c enddo print *,'icheckgrad=',icheckgrad goto (10,20,30) icheckgrad 10 call check_ecartint + call check_ecartint return 20 call check_cartgrad return diff --git a/source/unres/src_MD/COMMON.TORSION b/source/unres/src_MD/COMMON.TORSION index 6b6605f..6622cfa 100644 --- a/source/unres/src_MD/COMMON.TORSION +++ b/source/unres/src_MD/COMMON.TORSION @@ -9,11 +9,15 @@ C Torsional constants of the rotation about virtual-bond dihedral angles C 6/23/01 - constants for double torsionals double precision v1c,v1s,v2c,v2s integer ntermd_1,ntermd_2 - common /torsiond/ v1c(2,maxtermd_1,maxtor,maxtor,maxtor), - & v1s(2,maxtermd_1,maxtor,maxtor,maxtor), - & v2c(maxtermd_2,maxtermd_2,maxtor,maxtor,maxtor), - & v2s(maxtermd_2,maxtermd_2,maxtor,maxtor,maxtor), - & ntermd_1(maxtor,maxtor,maxtor),ntermd_2(maxtor,maxtor,maxtor) + common /torsiond/ + &v1c(2,maxtermd_1,maxtor,maxtor,maxtor), + &v1s(2,maxtermd_1,maxtor,maxtor,maxtor), + &v2c(maxtermd_2,maxtermd_2,maxtor,maxtor, + & maxtor), + &v2s(maxtermd_2,maxtermd_2,maxtor,maxtor, + & maxtor), + & ntermd_1(maxtor,maxtor,maxtor), + & ntermd_2(maxtor,maxtor,maxtor) C 9/18/99 - added Fourier coeffficients of the expansion of local energy C surface double precision b1,b2,cc,dd,ee,ctilde,dtilde,b2tilde diff --git a/source/unres/src_MD/parmread.F b/source/unres/src_MD/parmread.F index b7953b3..bfb4c22 100644 --- a/source/unres/src_MD/parmread.F +++ b/source/unres/src_MD/parmread.F @@ -210,6 +210,7 @@ C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203 C read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2, & ntheterm3,nsingle,ndouble +C print *, "tu" nntheterm=max0(ntheterm,ntheterm2,ntheterm3) read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1) do i=-ntyp1,-1 @@ -249,6 +250,7 @@ c VAR:ntethtyp is type of theta potentials type currently 0=glycine c VAR:1=non-glicyne non-proline 2=proline c VAR:negative values for D-aminoacid do i=0,nthetyp +C print *,i do j=-nthetyp,nthetyp do k=-nthetyp,nthetyp read (ithep,'(6a)',end=111,err=111) res1 @@ -642,15 +644,15 @@ c write (iout,*) 'ntortyp',ntortyp C C 6/23/01 Read parameters for double torsionals C - do i=0,ntortyp-1 - do j=-ntortyp+1,ntortyp-1 - do k=-ntortyp+1,ntortyp-1 + do i=1,ntortyp + do j=1,ntortyp + do k=1,ntortyp read (itordp,'(3a1)',end=114,err=114) t1,t2,t3 -c write (iout,*) "OK onelett", -c & i,j,k,t1,t2,t3 + write (iout,*) "OK onelett", + & i,j,k,t1,t2,t3 - if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) - & .or. t3.ne.toronelet(k)) then + if (t1.ne.onelett(i) .or. t2.ne.onelett(j) + & .or. t3.ne.onelett(k)) then write (iout,*) "Error in double torsional parameter file", & i,j,k,t1,t2,t3 #ifdef MPI @@ -729,6 +731,7 @@ cc maxinter is maximum interaction sites do j=1,nsccortyp read (isccor,*,end=1113,err=1113) nterm_sccor(i,j), & nlor_sccor(i,j) + print *,i,j,l v0ijsccor=0.0d0 v0ijsccor1=0.0d0 v0ijsccor2=0.0d0 @@ -772,8 +775,8 @@ cc maxinter is maximum interaction sites endif endif endif - read (isccor,*,end=1113,err=1113) kk,v1sccor(k,l,i,j) - & ,v2sccor(k,l,i,j) +C read (isccor,*,end=1113,err=1113) kk,v1sccor(k,l,i,j) +C & ,v2sccor(k,l,i,j) v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j) v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j) v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j) @@ -866,8 +869,8 @@ c b1(1,i)=0.0d0 c b1(2,i)=0.0d0 B1tilde(1,i) = b(3) B1tilde(2,i) =-b(5) - B1tilde(1,-i) =-b(3) - B1tilde(2,-i) =b(5) +C B1tilde(1,-i) =-b(3) +C B1tilde(2,-i) =b(5) c b1tilde(1,i)=0.0d0 c b1tilde(2,i)=0.0d0 B2(1,i) = b(2) diff --git a/source/unres/src_MD/ssMD.F b/source/unres/src_MD/ssMD.F index 9c4bfd3..6c7d523 100644 --- a/source/unres/src_MD/ssMD.F +++ b/source/unres/src_MD/ssMD.F @@ -589,7 +589,7 @@ cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss) & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR) call MPI_Gather(newnss,1,MPI_INTEGER, & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR) - displ(0)=0 +C displ(0)=0 do i=1,nfgtasks-1,1 displ(i)=i_newnss(i-1)+displ(i-1) enddo -- 1.7.9.5