intoduction of quartic restrains in multichain, bugfix in single chain

diff --git a/PARAM/bond_AM1_ext.parm b/PARAM/bond_AM1_ext.parm
index 333a43f..0f7bb5e 100644 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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