corrections to src_MD-M for the same energy as src_MD

diff --git a/source/unres/src_MD-M/COMMON.IOUNITS b/source/unres/src_MD-M/COMMON.IOUNITS
index a9ace0b..49b6db3 100644 (file)
--- a/source/unres/src_MD-M/COMMON.IOUNITS
+++ b/source/unres/src_MD-M/COMMON.IOUNITS
@@ -11,11 +11,11 @@ C General I/O units & files
       integer inp,iout,igeom,intin,ipdb,imol2,ipdbin,ithep,irotam,
      &        itorp,itordp,ifourier,ielep,isidep,iscpp,icbase,istat,
      &        ientin,ientout,izs1,isecpred,ibond,irest2,iifrag,icart,
-     &        irest1,isccor
+     &        irest1,isccor,ithep_pdb,irotam_pdb
       common /iounits/ inp,iout,igeom,intin,ipdb,imol2,ipdbin,ithep,
      &        irotam,itorp,itordp,ifourier,ielep,isidep,iscpp,icbase,
      &        istat,ientin,ientout,izs1,isecpred,ibond,irest2,iifrag,
-     &        icart,irest1,isccor
+     &        icart,irest1,isccor,ithep_pdb,irotam_pdb
       character*256 outname,intname,pdbname,mol2name,statname,intinname,
      &        entname,prefix,secpred,rest2name,qname,cartname,tmpdir,
      &        mremd_rst_name,curdir,pref_orig
@@ -38,9 +38,11 @@ C CSA I/O units & files
      & icsa_bank_reminimized,icsa_native_int,icsa_in,icsa_pdb
 C Parameter files
       character*256 bondname,thetname,rotname,torname,tordname,
-     &       fouriername,elename,sidename,scpname,sccorname,patname
+     &       fouriername,elename,sidename,scpname,sccorname,patname,
+     &       thetname_pdb,rotname_pdb
       common /parfiles/ bondname,thetname,rotname,torname,tordname,
-     &       fouriername,elename,sidename,scpname,sccorname,patname
+     &       fouriername,elename,sidename,scpname,sccorname,patname,
+     &       thetname_pdb,rotname_pdb
       character*3 pot
 C-----------------------------------------------------------------------
 C INP    - main input file

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 eb1cbd4..9ea982a 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -4546,7 +4546,8 @@ C
       logical lprn /.false./, lprn1 /.false./
       etheta=0.0D0
       do i=ithet_start,ithet_end
-        if (itype(i-1).eq.21) cycle
+        if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+     &(itype(i).eq.ntyp1)) cycle
         dethetai=0.0d0
         dephii=0.0d0
         dephii1=0.0d0
@@ -4556,7 +4557,8 @@ C
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
         enddo
-        if (i.gt.3 .and. itype(i-2).ne.21) then
+C        if (i.gt.3) then
+        if (i.gt.3 .and. itype(imax0(i-3,1)).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
@@ -4570,13 +4572,13 @@ C
           enddo
         else
           phii=0.0d0
-          ityp1=nthetyp+1
+          ityp1=ithetyp(itype(i-2))
           do k=1,nsingle
             cosph1(k)=0.0d0
             sinph1(k)=0.0d0
           enddo 
         endif
-        if (i.lt.nres .and. itype(i).ne.21) 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
@@ -4591,7 +4593,7 @@ C
           enddo
         else
           phii1=0.0d0
-          ityp3=nthetyp+1
+          ityp3=ithetyp(itype(i))
           do k=1,nsingle
             cosph2(k)=0.0d0
             sinph2(k)=0.0d0
@@ -4701,6 +4703,8 @@ C
      &   i,theta(i)*rad2deg,phii*rad2deg,
      &   phii1*rad2deg,ethetai
         etheta=etheta+ethetai
+        if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+     &      'ebend',i,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

diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F
index ae884df..37692d4 100644 (file)
--- a/source/unres/src_MD-M/initialize_p.F
+++ b/source/unres/src_MD-M/initialize_p.F
@@ -80,7 +80,9 @@ C
       igeom=  8
       intin=  9
       ithep= 11
+      ithep_pdb=51
       irotam=12
+      irotam_pdb=52
       itorp= 13
       itordp= 23
       ielep= 14

diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F
index 64203a2..183da9d 100644 (file)
--- a/source/unres/src_MD-M/parmread.F
+++ b/source/unres/src_MD-M/parmread.F
@@ -276,6 +276,17 @@ C
       enddo
       call flush(iout)
       endif
+      write (2,*) "Start reading THETA_PDB"
+      do i=1,ntyp
+        read (ithep_pdb,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
+     &    (bthet(j,i),j=1,2)
+        read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
+       read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
+       read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
+       sigc0(i)=sigc0(i)**2
+      enddo
+      write (2,*) "End reading THETA_PDB"
+      close (ithep_pdb)
 #endif
       close(ithep)
 #ifdef CRYST_SC
@@ -363,6 +374,48 @@ C
          enddo  
        endif
       enddo
+C
+C Read the parameters of the probability distribution/energy expression
+C of the side chains.
+C
+      do i=1,ntyp
+       read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
+        if (i.eq.10) then
+          dsc_inv(i)=0.0D0
+        else
+          dsc_inv(i)=1.0D0/dsc(i)
+        endif
+       if (i.ne.10) then
+        do j=1,nlob(i)
+          do k=1,3
+            do l=1,3
+              blower(l,k,j)=0.0D0
+            enddo
+          enddo
+        enddo  
+       bsc(1,i)=0.0D0
+        read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
+     &    ((blower(k,l,1),l=1,k),k=1,3)
+       do j=2,nlob(i)
+         read (irotam_pdb,*,end=112,err=112) bsc(j,i)
+         read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
+     &                                 ((blower(k,l,j),l=1,k),k=1,3)
+        enddo
+       do j=1,nlob(i)
+         do k=1,3
+           do l=1,k
+             akl=0.0D0
+             do m=1,3
+               akl=akl+blower(k,m,j)*blower(l,m,j)
+              enddo
+             gaussc(k,l,j,i)=akl
+             gaussc(l,k,j,i)=akl
+            enddo
+          enddo 
+       enddo
+       endif
+      enddo
+      close (irotam_pdb)
 #endif
       close(irotam)
 
@@ -590,23 +643,23 @@ cc maxinter is maximum interaction sites
 c      write (iout,*) 'ntortyp',ntortyp
       maxinter=3
 cc maxinter is maximum interaction sites
-      do l=1,maxinter
+      do l=1,maxinter    
       do i=1,nsccortyp
-        do j=1,nsccortyp
+       do j=1,nsccortyp
           read (isccor,*,end=113,err=113)
      & nterm_sccor(i,j),nlor_sccor(i,j)
           v0ijsccor=0.0d0
           si=-1.0d0
-
-          do k=1,nterm_sccor(i,j)
+  
+         do k=1,nterm_sccor(i,j)
             read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
-     &    ,v2sccor(k,l,i,j)
+     &    ,v2sccor(k,l,i,j) 
             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
             si=-si
           enddo
-          do k=1,nlor_sccor(i,j)
+         do k=1,nlor_sccor(i,j)
             read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
-     &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
+     &        vlor2sccor(k,i,j),vlor3sccor(k,i,j) 
             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
      &(1+vlor3sccor(k,i,j)**2)
           enddo
@@ -615,12 +668,12 @@ cc maxinter is maximum interaction sites
       enddo
       enddo
       close (isccor)
-
+      
 #endif      
       if (lprint) then
-        write (iout,'(/a/)') 'Torsional constants:'
-        do i=1,nsccortyp
-          do j=1,nsccortyp
+       write (iout,'(/a/)') 'Torsional constants:'
+       do i=1,nsccortyp
+         do j=1,nsccortyp
             write (iout,*) 'ityp',i,' jtyp',j
             write (iout,*) 'Fourier constants'
             do k=1,nterm_sccor(i,j)
@@ -628,7 +681,7 @@ cc maxinter is maximum interaction sites
             enddo
             write (iout,*) 'Lorenz constants'
             do k=1,nlor_sccor(i,j)
-              write (iout,'(3(1pe15.5))')
+             write (iout,'(3(1pe15.5))') 
      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
             enddo
           enddo

diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F
index 0b08ab7..bd17218 100644 (file)
--- a/source/unres/src_MD-M/readrtns_CSA.F
+++ b/source/unres/src_MD-M/readrtns_CSA.F
@@ -1862,8 +1862,16 @@ C Get parameter filenames and open the parameter files.
       open (ibond,file=bondname,status='old',readonly,shared)
       call getenv_loc('THETPAR',thetname)
       open (ithep,file=thetname,status='old',readonly,shared)
+#ifndef CRYST_THETA
+      call getenv_loc('THETPARPDB',thetname_pdb)
+      open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
+#endif
       call getenv_loc('ROTPAR',rotname)
       open (irotam,file=rotname,status='old',readonly,shared)
+#ifndef CRYST_SC
+      call getenv_loc('ROTPARPDB',rotname_pdb)
+      open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
+#endif
       call getenv_loc('TORPAR',torname)
       open (itorp,file=torname,status='old',readonly,shared)
       call getenv_loc('TORDPAR',tordname)
@@ -1888,9 +1896,17 @@ c      print *,"Processor",myrank," opened file IBOND"
       call getenv_loc('THETPAR',thetname)
       open (ithep,file=thetname,status='old',action='read')
 c      print *,"Processor",myrank," opened file ITHEP" 
+#ifndef CRYST_THETA
+      call getenv_loc('THETPARPDB',thetname_pdb)
+      open (ithep_pdb,file=thetname_pdb,status='old',action='read')
+#endif
       call getenv_loc('ROTPAR',rotname)
       open (irotam,file=rotname,status='old',action='read')
 c      print *,"Processor",myrank," opened file IROTAM" 
+#ifndef CRYST_SC
+      call getenv_loc('ROTPARPDB',rotname_pdb)
+      open (irotam_pdb,file=rotname_pdb,status='old',action='read')
+#endif
       call getenv_loc('TORPAR',torname)
       open (itorp,file=torname,status='old',action='read')
 c      print *,"Processor",myrank," opened file ITORP" 
@@ -1919,8 +1935,16 @@ C Get parameter filenames and open the parameter files.
       open (ibond,file=bondname,status='old')
       call getenv_loc('THETPAR',thetname)
       open (ithep,file=thetname,status='old')
+#ifndef CRYST_THETA
+      call getenv_loc('THETPARPDB',thetname_pdb)
+      open (ithep_pdb,file=thetname_pdb,status='old')
+#endif
       call getenv_loc('ROTPAR',rotname)
       open (irotam,file=rotname,status='old')
+#ifndef CRYST_SC
+      call getenv_loc('ROTPARPDB',rotname_pdb)
+      open (irotam_pdb,file=rotname_pdb,status='old')
+#endif
       call getenv_loc('TORPAR',torname)
       open (itorp,file=torname,status='old')
       call getenv_loc('TORDPAR',tordname)
@@ -1943,8 +1967,18 @@ C Get parameter filenames and open the parameter files.
       open (ibond,file=bondname,status='old',readonly)
       call getenv_loc('THETPAR',thetname)
       open (ithep,file=thetname,status='old',readonly)
+#ifndef CRYST_THETA
+      call getenv_loc('THETPARPDB',thetname_pdb)
+      print *,"thetname_pdb ",thetname_pdb
+      open (ithep_pdb,file=thetname_pdb,status='old',action='read')
+      print *,ithep_pdb," opened"
+#endif
       call getenv_loc('ROTPAR',rotname)
       open (irotam,file=rotname,status='old',readonly)
+#ifndef CRYST_SC
+      call getenv_loc('ROTPARPDB',rotname_pdb)
+      open (irotam_pdb,file=rotname_pdb,status='old',action='read')
+#endif
       call getenv_loc('TORPAR',torname)
       open (itorp,file=torname,status='old',readonly)
       call getenv_loc('TORDPAR',tordname)
@@ -2388,6 +2422,9 @@ C Initialize random number generator
 C
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
+#ifdef AMD64
+      integer*8 iseedi8
+#endif
 #ifdef MPI
       include 'mpif.h'
       logical OKRandom, prng_restart
@@ -2423,9 +2460,11 @@ C
       if (fg_rank.eq.0) then
       seed=seed*(me+1)+1
 #ifdef AMD64
-      if(me.eq.king)
-     &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
-      OKRandom = prng_restart(me,iseed)
+      iseedi8=dint(seed)
+      if(me.eq.king .or. .not. out1file)
+     &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
+      write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
+      OKRandom = prng_restart(me,iseedi8)
 #else
       do i=1,4
        tmp=65536.0d0**(4-i)