working cluster for nano parameters
[unres.git] / source / unres / src_MD-M / parmread.F
index e156c97..068897d 100644 (file)
@@ -101,6 +101,8 @@ c
         enddo
       endif
 C reading lipid parameters
+      write (iout,*) "iliptranpar",iliptranpar
+      call flush(iout)
        read(iliptranpar,*) pepliptran
        do i=1,ntyp
        read(iliptranpar,*) liptranene(i)
@@ -219,6 +221,8 @@ C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 C
       read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
      &  ntheterm3,nsingle,ndouble
+      write (iout,*) "ithep",ithep
+      call flush(iout)
       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
       read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
       do i=-ntyp1,-1
@@ -361,11 +365,12 @@ C Control printout of the coefficients of virtual-bond-angle potentials
 C
       if (lprint) then
         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
-        do i=1,nthetyp+1
-          do j=1,nthetyp+1
-            do k=1,nthetyp+1
+        do iblock=1,2
+        do i=0,nthetyp
+          do j=-nthetyp,nthetyp
+            do k=-nthetyp,nthetyp
               write (iout,'(//4a)') 
-     &         'Type ',onelett(i),onelett(j),onelett(k) 
+     &         'Type ',toronelet(i),toronelet(j),toronelet(k) 
               write (iout,'(//a,10x,a)') " l","a[l]"
               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
               write (iout,'(i2,1pe15.5)')
@@ -394,6 +399,7 @@ C
             enddo
           enddo
         enddo
+        enddo
       enddo
       call flush(iout)
       endif
@@ -624,6 +630,7 @@ C Read torsional parameters
 C
       read (itorp,*,end=113,err=113) ntortyp
       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
+      itortyp(ntyp1)=ntortyp
       do iblock=1,2
       do i=-ntyp,-1
        itortyp(i)=-itortyp(-i)
@@ -661,8 +668,9 @@ c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
       close (itorp)
       if (lprint) then
         write (iout,'(/a/)') 'Torsional constants:'
-        do i=1,ntortyp
-          do j=1,ntortyp
+        do iblock=1,2
+        do i=0,ntortyp-1
+          do j=-ntortyp+1,ntortyp-1
             write (iout,*) 'ityp',i,' jtyp',j
             write (iout,*) 'Fourier constants'
             do k=1,nterm(i,j,iblock)
@@ -676,6 +684,7 @@ c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
             enddo
           enddo
         enddo
+        enddo
       endif
 
 C
@@ -776,30 +785,79 @@ C Martix of D parameters for two dimesional fourier series
       endif
 #endif
 C read Czybyshev torsional parameters
-        read (itorkcc,*,end=121,err=121) nkcctyp
-        read (itorkcc,*,end=121,err=121) (itortyp_kcc(i),i=1,ntyp)
-        do i=-ntyp,-1
+      read (itorkcc,*,end=121,err=121) nkcctyp
+      read (itorkcc,*,end=121,err=121) (itortyp_kcc(i),i=1,ntyp)
+      do i=-ntyp,-1
         itortyp_kcc(i)=-itortyp_kcc(-i)
-        enddo
-        do i=1,nkcctyp
-         do j=1,nkcctyp
+      enddo
+      do i=0,nkcctyp
+        do j=0,nkcctyp
 C first we read the cos and sin gamma parameters
-          read (itorkcc,*,end=121,err=121) nterm_kcc(j,i)
-           do k=1,nterm_kcc(j,i)
-             read (itorkcc,*,end=121,err=121) 
-     &        v1_kcc(k,j,i),v2_kcc(k,j,i)
-           enddo
-C now Czybyshev parameters
-           read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
-           do k=1,nterm_kcc_Tb(j,i)
-             read (itorkcc,*,end=121,err=121) 
-     &        v1_chyb(k,j,i),v2_chyb(k,j,i)
-           enddo
+          read (itorkcc,*,end=121,err=121) 
+     &    nterm_kcc(j,i),nterm_kcc_Tb(j,i)
+C           read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
+          do k=1,nterm_kcc(j,i)
+            do l=1,nterm_kcc_Tb(j,i)
+              read (itorkcc,*,end=121,err=121) v11_chyb(l,k,j,i)
+            enddo
+            do l=1,nterm_kcc_Tb(j,i)
+              read (itorkcc,*,end=121,err=121) v21_chyb(l,k,j,i)
+            enddo
+            do l=1,nterm_kcc_Tb(j,i)
+              read (itorkcc,*,end=121,err=121) v12_chyb(l,k,j,i)
+            enddo
+            do l=1,nterm_kcc_Tb(j,i)
+              read (itorkcc,*,end=121,err=121) v22_chyb(l,k,j,i)
+            enddo
+            read (itorkcc,*,end=121,err=121) v1_kcc(k,j,i)
+            read (itorkcc,*,end=121,err=121) v2_kcc(k,j,i)
           enddo
-         enddo
+        enddo
+      enddo
+      if (lprint) then
+c Print valence-torsional parameters
+        write (iout,'(a)') 
+     &    "Parameters of the valence-torsional potentials"
+        do i=0,nkcctyp
+        do j=0,nkcctyp
+        write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
+        write (iout,'(2a20,a15)') "v_kcc","v1_chyb","v2_chyb"
+        do k=1,nterm_kcc(j,i)
+          write (iout,'(i5,f15.10,i5,2f15.10)') 
+     &      k,v1_kcc(k,j,i),1,v11_chyb(1,k,j,i),v21_chyb(1,k,j,i)
+          do l=2,nterm_kcc_Tb(j,i)
+            write (iout,'(20x,i5,2f15.10)') 
+     &        l,v11_chyb(l,k,j,i),v21_chyb(l,k,j,i)
+          enddo
+          write (iout,'(i5,f15.10,i5,2f15.10)') 
+     &      k,v2_kcc(k,j,i),1,v12_chyb(1,k,j,i),v22_chyb(1,k,j,i)
+          do l=2,nterm_kcc_Tb(j,i)
+            write (iout,'(20x,i5,2f15.10)') 
+     &        l,v12_chyb(l,k,j,i),v22_chyb(l,k,j,i)
+          enddo
+          write (iout,'(a)')
+        enddo
+        enddo
+        enddo
+      endif
 C here will be the apropriate recalibrating for D-aminoacid
-
-
+C        read (ithetkcc,*,end=121,err=121) nkcctyp
+      do i=0,nkcctyp
+        read (ithetkcc,*,end=121,err=121) nbend_kcc_Tb(i)
+        do j=1,nbend_kcc_Tb(i)
+          read (ithetkcc,*,end=121,err=121) v1bend_chyb(j,i)
+        enddo
+      enddo
+      if (lprint) then
+        write (iout,'(a)') 
+     &    "Parameters of the valence-only potentials"
+        do i=0,nkcctyp
+        write (iout,'(2a)') "Type ",toronelet(i)
+        do k=1,nbend_kcc_Tb(i)
+          write(iout,'(i5,f15.10)') k,v1bend_chyb(k,i)
+        enddo
+        enddo
+      endif
 C Read of Side-chain backbone correlation parameters
 C Modified 11 May 2012 by Adasko
 CCC
@@ -885,13 +943,48 @@ cc maxinter is maximum interaction sites
 #else
       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
 c      write (iout,*) 'ntortyp',ntortyp
+      do i=-ntyp,-1
+        isccortyp(i)=-isccortyp(-i)
+      enddo
       maxinter=3
 cc maxinter is maximum interaction sites
+C NOW READING FOR LL aminoacid chirality
       do l=1,maxinter
       do i=1,nsccortyp
         do j=1,nsccortyp
           read (isccor,*,end=119,err=119)
      & nterm_sccor(i,j),nlor_sccor(i,j)
+          nterm_sccor(-i,-j)=nterm_sccor(i,j)
+          v0ijsccor=0.0d0
+          si=-1.0d0
+
+          do k=1,nterm_sccor(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
+            v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
+            v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
+
+          enddo
+          do k=1,nlor_sccor(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)
+          enddo
+          v0sccor(l,i,j)=v0ijsccor
+          v0sccor(l,-i,-j)=v0ijsccor
+        enddo
+      enddo
+      enddo
+C NOW READING FOR LD aminoacid chirality
+      do l=1,maxinter
+      do i=1,nsccortyp
+        do j=-nsccortyp,-1
+          read (isccor,*,end=119,err=119)
+     & nterm_sccor(i,j),nlor_sccor(i,j)
+          nterm_sccor(-i,-j)=nterm_sccor(i,j)
           v0ijsccor=0.0d0
           si=-1.0d0
 
@@ -900,6 +993,9 @@ cc maxinter is maximum interaction sites
      &    ,v2sccor(k,l,i,j)
             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
             si=-si
+            v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
+            v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
+
           enddo
           do k=1,nlor_sccor(i,j)
             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
@@ -907,15 +1003,30 @@ cc maxinter is maximum interaction sites
             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
      &(1+vlor3sccor(k,i,j)**2)
           enddo
-          v0sccor(i,j,iblock)=v0ijsccor
+          v0sccor(l,i,j)=v0ijsccor
+          v0sccor(l,-i,-j)=v0ijsccor
         enddo
       enddo
       enddo
+C symetry for GLY
+             i=10
+             do l=1,maxinter
+                do j=-nsccortyp,-1
+                do k=1,nterm_sccor(i,j)
+                v1sccor(k,l,10,j)=v1sccor(k,l,10,-j)
+                v2sccor(k,l,10,j)=-v2sccor(k,l,10,-j)
+                v1sccor(k,l,j,10)=v1sccor(k,l,-j,10)
+                v2sccor(k,l,j,10)=-v2sccor(k,l,-j,10)
+               enddo
+              enddo
+             enddo
+
       close (isccor)
 
 #endif      
       if (lprint) then
         write (iout,'(/a/)') 'Torsional constants:'
+        do l=1,maxinter
         do i=1,nsccortyp
           do j=1,nsccortyp
             write (iout,*) 'ityp',i,' jtyp',j
@@ -930,6 +1041,7 @@ cc maxinter is maximum interaction sites
             enddo
           enddo
         enddo
+        enddo
       endif
 
 C
@@ -941,7 +1053,29 @@ C
         write (iout,*) "Coefficients of the cumulants"
       endif
       read (ifourier,*) nloctyp
-
+#ifdef NEWCORR
+      read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
+      read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
+      itype2loc(ntyp1)=nloctyp
+      iloctyp(nloctyp)=ntyp1
+#else
+      do i=1,ntyp1
+        itype2loc(i)=itortyp(i)
+      enddo
+      iloctyp(0)=10
+      iloctyp(1)=9
+      iloctyp(2)=20
+      iloctyp(3)=ntyp1
+#endif
+      do i=1,ntyp1
+        itype2loc(-i)=-itype2loc(i)
+      enddo
+      do i=1,nloctyp
+        iloctyp(-i)=-iloctyp(i)
+      enddo
+      write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
+      write (iout,*) "nloctyp",nloctyp,
+     &  " iloctyp",(iloctyp(i),i=0,nloctyp)
       do i=0,nloctyp-1
         read (ifourier,*,end=115,err=115)
         read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
@@ -956,30 +1090,34 @@ C
         write (iout,*) 'Type',i
         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
         endif
-c        B1(1,i)  = b(3)
-c        B1(2,i)  = b(5)
-c        B1(1,-i) = b(3)
-c        B1(2,-i) = -b(5)
+        b(3,-i)= b(3,i)
+        b(5,-i)=-b(5,i)
+        b(4,-i)=-b(4,i)
+        b(2,-i)= b(2,i)     
+        B1(1,i)  = b(3,i)
+        B1(2,i)  = b(5,i)
+        B1(1,-i) = b(3,i)
+        B1(2,-i) = -b(5,i)
 c        b1(1,i)=0.0d0
 c        b1(2,i)=0.0d0
-c        B1tilde(1,i) = b(3)
-c        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
-c        B2(1,i)  = b(2)
-c        B2(2,i)  = b(4)
-c        B2(1,-i)  =b(2)
-c        B2(2,-i)  =-b(4)
         B1tilde(1,i) = b(3,i)
         B1tilde(2,i) =-b(5,i)
-C        B1tilde(1,-i) =-b(3,i)
-C        B1tilde(2,-i) =b(5,i)
-        b1tilde(1,i)=0.0d0
-        b1tilde(2,i)=0.0d0
+        B1tilde(1,-i) =-b(3,i)
+        B1tilde(2,-i) =b(5,i)
+c        b1tilde(1,i)=0.0d0
+c        b1tilde(2,i)=0.0d0
         B2(1,i)  = b(2,i)
         B2(2,i)  = b(4,i)
+        B2(1,-i)  =b(2,i)
+        B2(2,-i)  =-b(4,i)
+cc        B1tilde(1,i) = b(3,i)
+cc        B1tilde(2,i) =-b(5,i)
+C        B1tilde(1,-i) =-b(3,i)
+C        B1tilde(2,-i) =b(5,i)
+cc        b1tilde(1,i)=0.0d0
+cc        b1tilde(2,i)=0.0d0
+cc        B2(1,i)  = b(2,i)
+cc        B2(2,i)  = b(4,i)
 C        B2(1,-i)  =b(2,i)
 C        B2(2,-i)  =-b(4,i)
 
@@ -1051,7 +1189,7 @@ c        ee(2,1,i)=0.0d0
 c        ee(1,2,i)=0.0d0
 c        ee(2,1,i)=ee(1,2,i)
       enddo
-c      lprint=.true.
+      lprint=.true.
       if (lprint) then
       do i=1,nloctyp
         write (iout,*) 'Type',i
@@ -1073,7 +1211,7 @@ c      lprint=.true.
         enddo
       enddo
       endif
-c      lprint=.false.
+      lprint=.false.
 
 C 
 C Read electrostatic-interaction parameters
@@ -1155,6 +1293,8 @@ C---------------------- GB or BP potential -----------------------------
 C now we start reading lipid
       do i=1,ntyp
        read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
+       write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
+
 C       print *,"WARNING!!"
 C       do j=1,ntyp
 C       epslip(i,j)=epslip(i,j)+0.05d0
@@ -1276,6 +1416,26 @@ c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
          endif
         enddo
       enddo
+      write(iout,*) "tube param"
+      read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep,
+     & ccavtubpep,dcavtubpep,tubetranenepep
+      sigmapeptube=sigmapeptube**6
+      sigeps=dsign(1.0D0,epspeptube)
+      epspeptube=dabs(epspeptube)
+      pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
+      pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
+      write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
+      do i=1,ntyp
+       read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i),
+     & ccavtub(i),dcavtub(i),tubetranene(i)
+       sigmasctube=sigmasctube**6
+       sigeps=dsign(1.0D0,epssctube)
+       epssctube=dabs(epssctube)
+       sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
+       sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
+      write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
+      enddo
+
 #ifdef OLDSCP
 C
 C Define the SC-p interaction constants (hard-coded; old style)