Clean-up in parmread for src-MD and src-MD-M
[unres.git] / source / unres / src_MD-M / parmread.F
index 44998e3..154fcf4 100644 (file)
@@ -243,12 +243,21 @@ C
         enddo
       enddo
       enddo
+c VAR:iblock means terminally blocking group 1=non-proline 2=proline
       do iblock=1,2 
+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
         do j=-nthetyp,nthetyp
           do k=-nthetyp,nthetyp
             read (ithep,'(6a)',end=111,err=111) res1
             read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
+c VAR: aa0thet is variable describing the average value of Foureir
+c VAR: expansion series
+c VAR: aathet is foureir expansion in theta/2 angle for full formula
+c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
+Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
             read (ithep,*,end=111,err=111)
      &(aathet(l,i,j,k,iblock),l=1,ntheterm)
             read (ithep,*,end=111,err=111)
@@ -269,7 +278,8 @@ C
 C
 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 C coefficients of theta-and-gamma-dependent terms are zero.
-C
+C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
+C RECOMENTDED AFTER VERSION 3.3)
 c      do i=1,nthetyp
 c        do j=1,nthetyp
 c          do l=1,ntheterm
@@ -285,6 +295,7 @@ c        enddo
 c        aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
 c      enddo
 c      enddo
+C AND COMMENT THE LOOPS BELOW
       do i=1,nthetyp
         do j=1,nthetyp
           do l=1,ntheterm
@@ -300,7 +311,7 @@ c      enddo
         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
       enddo
       enddo
-
+C TILL HERE
 C Substitution for D aminoacids from symmetry.
       do iblock=1,2
       do i=-nthetyp,0
@@ -377,6 +388,51 @@ C
       enddo
       call flush(iout)
       endif
+      write (2,*) "Start reading THETA_PDB"
+      do i=1,ntyp
+        read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
+     &    (bthet(j,i,1,1),j=1,2)
+        read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
+        read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
+        read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
+        sigc0(i)=sigc0(i)**2
+      enddo
+      do i=1,ntyp
+      athet(1,i,1,-1)=athet(1,i,1,1)
+      athet(2,i,1,-1)=athet(2,i,1,1)
+      bthet(1,i,1,-1)=-bthet(1,i,1,1)
+      bthet(2,i,1,-1)=-bthet(2,i,1,1)
+      athet(1,i,-1,1)=-athet(1,i,1,1)
+      athet(2,i,-1,1)=-athet(2,i,1,1)
+      bthet(1,i,-1,1)=bthet(1,i,1,1)
+      bthet(2,i,-1,1)=bthet(2,i,1,1)
+      enddo
+      do i=-ntyp,-1
+      a0thet(i)=a0thet(-i)
+      athet(1,i,-1,-1)=athet(1,-i,1,1)
+      athet(2,i,-1,-1)=-athet(2,-i,1,1)
+      bthet(1,i,-1,-1)=bthet(1,-i,1,1)
+      bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
+      athet(1,i,-1,1)=athet(1,-i,1,1)
+      athet(2,i,-1,1)=-athet(2,-i,1,1)
+      bthet(1,i,-1,1)=-bthet(1,-i,1,1)
+      bthet(2,i,-1,1)=bthet(2,-i,1,1)
+      athet(1,i,1,-1)=-athet(1,-i,1,1)
+      athet(2,i,1,-1)=athet(2,-i,1,1)
+      bthet(1,i,1,-1)=bthet(1,-i,1,1)
+      bthet(2,i,1,-1)=-bthet(2,-i,1,1)
+      theta0(i)=theta0(-i)
+      sig0(i)=sig0(-i)
+      sigc0(i)=sigc0(-i)
+       do j=0,3
+        polthet(j,i)=polthet(j,-i)
+       enddo
+       do j=1,3
+         gthet(j,i)=gthet(j,-i)
+       enddo
+      enddo
+      write (2,*) "End reading THETA_PDB"
+      close (ithep_pdb)
 #endif
       close(ithep)
 #ifdef CRYST_SC
@@ -470,7 +526,6 @@ C Read scrot parameters for potentials determined from all-atom AM1 calculations
 C added by Urszula Kozlowska 07/11/2007
 C
       do i=1,ntyp
-CC    TU JEST ZLE musibyc ntyp
        read (irotam,*,end=112,err=112) 
        if (i.eq.10) then 
          read (irotam,*,end=112,err=112) 
@@ -480,6 +535,48 @@ CC    TU JEST ZLE musibyc ntyp
          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)
 
@@ -734,7 +831,7 @@ cc maxinter is maximum interaction sites
             si=-si
           enddo
           do k=1,nlor_sccor(i,j)
-            read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
+            read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
      &(1+vlor3sccor(k,i,j)**2)
@@ -748,26 +845,26 @@ cc maxinter is maximum interaction sites
       enddo
       close (isccor)
 #else
-      read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
+      read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
 c      write (iout,*) 'ntortyp',ntortyp
       maxinter=3
 cc maxinter is maximum interaction sites
       do l=1,maxinter
       do i=1,nsccortyp
         do j=1,nsccortyp
-          read (isccor,*,end=113,err=113)
+          read (isccor,*,end=119,err=119)
      & nterm_sccor(i,j),nlor_sccor(i,j)
           v0ijsccor=0.0d0
           si=-1.0d0
 
           do k=1,nterm_sccor(i,j)
-            read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
+            read (isccor,*,end=119,err=119) kk,v1sccor(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)
-            read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
+            read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
      &(1+vlor3sccor(k,i,j)**2)