Dzialajaczy dla D aminokwasow whama(mono) psuedosymetria dla GABowego
[unres.git] / source / unres / src_MD / parmread.F
index 654c305..9d5100b 100644 (file)
@@ -28,6 +28,7 @@ C
       include 'COMMON.SETUP'
       character*1 t1,t2,t3
       character*1 onelett(4) /"G","A","P","D"/
+      character*1 toronelet(-2:2) /"p","a","G","A","P"/
       logical lprint,LaTeX
       dimension blower(3,3,maxlob)
       dimension b(13)
@@ -102,13 +103,47 @@ C Read the parameters of the probability distribution/energy expression
 C of the virtual-bond valence angles theta
 C
       do i=1,ntyp
-        read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
-     &    (bthet(j,i),j=1,2)
+        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
       close (ithep)
       if (lprint) then
       if (.not.LaTeX) then
@@ -119,7 +154,7 @@ C
      & '        B1    ','         B2   '        
         do i=1,ntyp
           write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
-     &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
+     &        a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
         enddo
         write (iout,'(/a/9x,5a/79(1h-))') 
      & 'Parameters of the expression for sigma(theta_c):',
@@ -146,7 +181,8 @@ C
      & '   b1*10^1    ','    b2*10^1   '        
         do i=1,ntyp
           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
-     &        a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
+     &        a0thet(i),(100*athet(j,i,1,1),j=1,2),
+     &        (10*bthet(j,i,1,1),j=1,2)
         enddo
        write (iout,'(/a/9x,5a/79(1h-))') 
      & 'Parameters of the expression for sigma(theta_c):',
@@ -472,7 +508,7 @@ C
       enddo
 c      write (iout,*) 'ntortyp',ntortyp
       do i=0,ntortyp-1
-       do j=-ntortyp,ntortyp
+       do j=-ntortyp+1,ntortyp-1
          read (itorp,*,end=113,err=113) nterm(i,j,iblock),
      &          nlor(i,j,iblock)
           nterm(-i,-j,iblock)=nterm(i,j,iblock)
@@ -484,8 +520,11 @@ c      write (iout,*) 'ntortyp',ntortyp
      &      v2(k,i,j,iblock) 
             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
-            v0ij=v0ij+si*v1(k,i,j)
+            v0ij=v0ij+si*v1(k,i,j,iblock)
             si=-si
+c         write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
+c         write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
+c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
           enddo
          do k=1,nlor(i,j,iblock)
             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
@@ -504,11 +543,12 @@ c      write (iout,*) 'ntortyp',ntortyp
          do j=1,ntortyp
             write (iout,*) 'ityp',i,' jtyp',j
             write (iout,*) 'Fourier constants'
-            do k=1,nterm(i,j)
-             write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
+            do k=1,nterm(i,j,iblock)
+             write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
+     &        v2(k,i,j,iblock)
             enddo
             write (iout,*) 'Lorenz constants'
-            do k=1,nlor(i,j)
+            do k=1,nlor(i,j,iblock)
              write (iout,'(3(1pe15.5))') 
      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
             enddo
@@ -518,12 +558,16 @@ c      write (iout,*) 'ntortyp',ntortyp
 C
 C 6/23/01 Read parameters for double torsionals
 C
-      do i=1,ntortyp
-        do j=1,ntortyp
-          do k=1,ntortyp
+      do iblock=1,2
+      do i=0,ntortyp-1
+        do j=-ntortyp+1,ntortyp-1
+          do k=-ntortyp+1,ntortyp-1
             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
-            if (t1.ne.onelett(i) .or. t2.ne.onelett(j) 
-     &        .or. t3.ne.onelett(k)) then
+c              write (iout,*) "OK onelett",
+c     &         i,j,k,t1,t2,t3
+
+            if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) 
+     &        .or. t3.ne.toronelet(k)) then
               write (iout,*) "Error in double torsional parameter file",
      &         i,j,k,t1,t2,t3
 #ifdef MPI
@@ -531,54 +575,81 @@ C
 #endif
                stop "Error in double torsional parameter file"
             endif
-            read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
-     &         ntermd_2(i,j,k)
-            read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
-     &         ntermd_1(i,j,k))
-            read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
-     &         ntermd_1(i,j,k))
-            read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
-     &         ntermd_1(i,j,k))
-            read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
-     &         ntermd_1(i,j,k))
-            read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
-     &         v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
-     &         m=1,l-1),l=1,ntermd_2(i,j,k))
-          enddo
-        enddo
-      enddo
+            read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
+     &         ntermd_2(i,j,k,iblock)
+            ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
+            ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
+            read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
+     &         ntermd_1(i,j,k,iblock))
+            read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
+     &         ntermd_1(i,j,k,iblock))
+            read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
+     &         ntermd_1(i,j,k,iblock))
+            read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
+     &         ntermd_1(i,j,k,iblock))
+C Martix of D parameters for one dimesional foureir series
+            do l=1,ntermd_1(i,j,k,iblock)
+             v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
+             v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
+             v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
+             v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
+c            write(iout,*) "whcodze" ,
+c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
+            enddo
+            read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
+     &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
+     &         v2s(m,l,i,j,k,iblock),
+     &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
+C Martix of D parameters for two dimesional fourier series
+            do l=1,ntermd_2(i,j,k,iblock)
+             do m=1,l-1
+             v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
+             v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
+             v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
+             v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
+             enddo!m
+            enddo!l
+          enddo!k
+        enddo!j
+      enddo!i
+      enddo!iblock
       if (lprint) then
       write (iout,*) 
       write (iout,*) 'Constants for double torsionals'
+      do iblock=1,2
       do i=1,ntortyp
-        do j=1,ntortyp 
-          do k=1,ntortyp
+        do j=-ntortyp,ntortyp 
+          do k=-ntortyp,ntortyp
             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
-     &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
+     &        ' nsingle',ntermd_1(i,j,k,iblock),
+     &        ' ndouble',ntermd_2(i,j,k,iblock)
             write (iout,*)
             write (iout,*) 'Single angles:'
-            do l=1,ntermd_1(i,j,k)
-              write (iout,'(i5,2f10.5,5x,2f10.5)') l,
-     &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
-     &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
+            do l=1,ntermd_1(i,j,k,iblock)
+              write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
+     &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
+     &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
+     &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
             enddo
             write (iout,*)
             write (iout,*) 'Pairs of angles:'
-            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
-            do l=1,ntermd_2(i,j,k)
+            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
+            do l=1,ntermd_2(i,j,k,iblock)
               write (iout,'(i5,20f10.5)') 
-     &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
+     &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
             enddo
             write (iout,*)
-            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
-            do l=1,ntermd_2(i,j,k)
+            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
+            do l=1,ntermd_2(i,j,k,iblock)
               write (iout,'(i5,20f10.5)') 
-     &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
+     &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
+     &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
             enddo
             write (iout,*)
           enddo
         enddo
       enddo
+      enddo
       endif
 #endif
 C Read of Side-chain backbone correlation parameters
@@ -588,9 +659,10 @@ C
       read (isccor,*,end=113,err=113) nsccortyp
 #ifdef SCCORPDB
       read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
-#else
-      read (isccor,*,end=113,err=113) (isccortyp(i),i=-ntyp,ntyp)
-#endif
+      do i=-ntyp,-1
+        isccortyp(i)=-isccortyp(-i)
+      enddo
+      iscprol=isccortyp(20)
 c      write (iout,*) 'ntortyp',ntortyp
       maxinter=3
 cc maxinter is maximum interaction sites
@@ -603,7 +675,24 @@ cc maxinter is maximum interaction sites
   
          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)
+            if (j.eq.iscprol) then
+             v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
+     &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
+             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
+     &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
+             v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
+             v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
+             v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
+             v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)          
+            else
+             v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
+             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
+             v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
+             v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
+             v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
+             v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
+            endif 
             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
             si=-si
           enddo
@@ -618,7 +707,38 @@ cc maxinter is maximum interaction sites
       enddo
       enddo
       close (isccor)
-      
+#else
+      read (isccor,*,end=113,err=113) (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)
+      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)
+     &    ,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),
+     &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
+            v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
+     &(1+vlor3sccor(k,i,j)**2)
+          enddo
+          v0sccor(i,j)=v0ijsccor
+        enddo
+      enddo
+      enddo
+      close (isccor)
+
+#endif      
       if (lprint) then
        write (iout,'(/a/)') 'Torsional constants:'
        do i=1,nsccortyp
@@ -662,7 +782,7 @@ 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) 
+        B1tilde(2,-i) =b(5)
 c        b1tilde(1,i)=0.0d0
 c        b1tilde(2,i)=0.0d0
         B2(1,i)  = b(2)