unres_package_Oct_2016 from emilial
[unres4.git] / source / unres / io_config.f90
index 0b10e11..490ecad 100644 (file)
@@ -5,6 +5,7 @@
       use io_base
       use geometry_data
       use geometry
+      use control_data, only:maxterm_sccor
       implicit none
 !-----------------------------------------------------------------------------
 ! Max. number of residue types and parameters in expressions for 
@@ -21,7 +22,7 @@
 !      parameter (maxtor=4,maxterm=10)
 !-----------------------------------------------------------------------------
 ! Max number of torsional terms in SCCOR
-      integer,parameter :: maxterm_sccor=6
+!el      integer,parameter :: maxterm_sccor=6
 !-----------------------------------------------------------------------------
       character(len=1),dimension(:),allocatable :: secstruc    !(maxres)
 !-----------------------------------------------------------------------------
@@ -29,7 +30,7 @@
 !
 !-----------------------------------------------------------------------------
       contains
-#ifndef WHAM_RUN
+#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
 !-----------------------------------------------------------------------------
 ! bank.F    io_csa
 !-----------------------------------------------------------------------------
         endif        
       enddo
       ndih_constr=ii
-      deallocate(secstruc)
+!      deallocate(secstruc)
       return
 100   continue
       write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
-      deallocate(secstruc)
+!      deallocate(secstruc)
       return
 110   continue
       write(iout,'(A20)')'Error reading FTORS'
-      deallocate(secstruc)
+!      deallocate(secstruc)
       return
       end subroutine secstrp2dihc
 !-----------------------------------------------------------------------------
        ' in position',i4)
       return
       end subroutine read_secstr_pred
-#endif
+!#endif
 !-----------------------------------------------------------------------------
 ! parmread.F
 !-----------------------------------------------------------------------------
 !      real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
 !el      allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
 !el      allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
-
 !
 ! For printing parameters after they are read set the following in the UNRES
 ! C-shell script:
       allocate(restok(ntyp+1)) !(ntyp+1)
       allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 
+      dsc(:)=0.0d0
+      dsc_inv(:)=0.0d0
+
 #ifdef CRYST_BOND
       read (ibond,*) vbldp0,akp,mp,ip,pstok
       do i=1,ntyp
       do i=-ntyp,-1
        itortyp(i)=-itortyp(-i)
       enddo
-      write (iout,*) 'ntortyp',ntortyp
+!      itortyp(ntyp1)=ntortyp
+!      itortyp(-ntyp1)=-ntortyp
       do iblock=1,2 
+      write (iout,*) 'ntortyp',ntortyp
       do i=0,ntortyp-1
         do j=-ntortyp+1,ntortyp-1
           read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
 
             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
               .or. t3.ne.toronelet(k)) then
-              write (iout,*) "Error in double torsional parameter file",&
+             write (iout,*) "Error in double torsional parameter file",&
                i,j,k,t1,t2,t3
 #ifdef MPI
               call MPI_Finalize(Ierror)
         write (iout,*) "Coefficients of the cumulants"
       endif
       read (ifourier,*) nloctyp
-write(iout,*) "nloctyp",nloctyp
+!write(iout,*) "nloctyp",nloctyp
 !el from module energy-------
       allocate(b1(2,-nloctyp-1:nloctyp+1))     !(2,-maxtor:maxtor)
       allocate(b2(2,-nloctyp-1:nloctyp+1))     !(2,-maxtor:maxtor)
@@ -1807,14 +1812,17 @@ write(iout,*) "nloctyp",nloctyp
       allocate(ee(2,2,-nloctyp-1:nloctyp+1))
       allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
       allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
+! el
+        b1(1,:)=0.0d0
+        b1(2,:)=0.0d0
 !--------------------------------
 
       do i=0,nloctyp-1
         read (ifourier,*,end=115,err=115)
         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
         if (lprint) then
-        write (iout,*) 'Type',i
-        write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
+          write (iout,*) 'Type',i
+          write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
         endif
         B1(1,i)  = b(3)
         B1(2,i)  = b(5)
@@ -1921,7 +1929,6 @@ write(iout,*) "nloctyp",nloctyp
         enddo
       enddo
       endif
-
 ! 
 ! Read electrostatic-interaction parameters
 !
@@ -1984,72 +1991,85 @@ write(iout,*) "nloctyp",nloctyp
       if(me.eq.king) &
        write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
        ', exponents are ',expon,2*expon 
-      goto (10,20,30,30,40) ipot
+!      goto (10,20,30,30,40) ipot
+      select case(ipot)
 !----------------------- LJ potential ---------------------------------
-   10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
-         (sigma0(i),i=1,ntyp)
-      if (lprint) then
-       write (iout,'(/a/)') 'Parameters of the LJ potential:'
-       write (iout,'(a/)') 'The epsilon array:'
-       call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
-       write (iout,'(/a)') 'One-body parameters:'
-       write (iout,'(a,4x,a)') 'residue','sigma'
-       write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
-      endif
-      goto 50
+       case (1)
+!   10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+         read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+           (sigma0(i),i=1,ntyp)
+        if (lprint) then
+         write (iout,'(/a/)') 'Parameters of the LJ potential:'
+         write (iout,'(a/)') 'The epsilon array:'
+         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
+         write (iout,'(/a)') 'One-body parameters:'
+         write (iout,'(a,4x,a)') 'residue','sigma'
+         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
+        endif
+!      goto 50
 !----------------------- LJK potential --------------------------------
-   20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
-        (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
-      if (lprint) then
-       write (iout,'(/a/)') 'Parameters of the LJK potential:'
-       write (iout,'(a/)') 'The epsilon array:'
-       call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
-       write (iout,'(/a)') 'One-body parameters:'
-       write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
-        write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
-              i=1,ntyp)
-      endif
-      goto 50
+       case(2)
+!   20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+         read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+          (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
+        if (lprint) then
+         write (iout,'(/a/)') 'Parameters of the LJK potential:'
+         write (iout,'(a/)') 'The epsilon array:'
+         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
+         write (iout,'(/a)') 'One-body parameters:'
+         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
+          write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
+                i=1,ntyp)
+        endif
+!      goto 50
 !---------------------- GB or BP potential -----------------------------
-   30 do i=1,ntyp
-       read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
-      enddo
-      read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
-      read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
-      read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
-      read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
-! For the GB potential convert sigma'**2 into chi'
-      if (ipot.eq.4) then
-       do i=1,ntyp
-         chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
+       case(3:4)
+!   30 do i=1,ntyp
+        do i=1,ntyp
+         read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
         enddo
-      endif
-      if (lprint) then
-       write (iout,'(/a/)') 'Parameters of the BP potential:'
-       write (iout,'(a/)') 'The epsilon array:'
-       call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
-       write (iout,'(/a)') 'One-body parameters:'
-       write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',&
-             '    chip  ','    alph  '
-       write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
-                           chip(i),alp(i),i=1,ntyp)
-      endif
-      goto 50
+        read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
+        read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
+        read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
+        read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
+! For the GB potential convert sigma'**2 into chi'
+        if (ipot.eq.4) then
+         do i=1,ntyp
+           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
+          enddo
+        endif
+        if (lprint) then
+         write (iout,'(/a/)') 'Parameters of the BP potential:'
+         write (iout,'(a/)') 'The epsilon array:'
+         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
+         write (iout,'(/a)') 'One-body parameters:'
+         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',&
+               '    chip  ','    alph  '
+         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
+                             chip(i),alp(i),i=1,ntyp)
+        endif
+!      goto 50
 !--------------------- GBV potential -----------------------------------
-   40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
-        (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
-        (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
-      if (lprint) then
-       write (iout,'(/a/)') 'Parameters of the GBV potential:'
-       write (iout,'(a/)') 'The epsilon array:'
-       call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
-       write (iout,'(/a)') 'One-body parameters:'
-       write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',&
-            's||/s_|_^2','    chip  ','    alph  '
-       write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
-                 sigii(i),chip(i),alp(i),i=1,ntyp)
-      endif
-   50 continue
+       case(5)
+!   40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+        read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+          (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
+          (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
+        if (lprint) then
+          write (iout,'(/a/)') 'Parameters of the GBV potential:'
+          write (iout,'(a/)') 'The epsilon array:'
+          call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
+          write (iout,'(/a)') 'One-body parameters:'
+          write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',&
+              's||/s_|_^2','    chip  ','    alph  '
+          write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
+                   sigii(i),chip(i),alp(i),i=1,ntyp)
+        endif
+       case default
+        write(iout,*)"Wrong ipot"
+!   50 continue
+      end select
+      continue
       close (isidep)
 !-----------------------------------------------------------------------
 ! Calculate the "working" parameters of SC interactions.
@@ -2058,19 +2078,19 @@ write(iout,*) "nloctyp",nloctyp
       allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
       allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
       do i=1,ntyp1
-       do j=1,ntyp1
-         aa(i,j)=0.0D0
-         bb(i,j)=0.0D0
-         chi(i,j)=0.0D0
-         sigma(i,j)=0.0D0
-         r0(i,j)=0.0D0
+        do j=1,ntyp1
+          aa(i,j)=0.0D0
+          bb(i,j)=0.0D0
+          chi(i,j)=0.0D0
+          sigma(i,j)=0.0D0
+          r0(i,j)=0.0D0
         enddo
       enddo
 !--------------------------------
 
       do i=2,ntyp
         do j=1,i-1
-         eps(i,j)=eps(j,i)
+          eps(i,j)=eps(j,i)
         enddo
       enddo
       do i=1,ntyp
@@ -2259,88 +2279,7 @@ write(iout,*) "nloctyp",nloctyp
       stop
       return
       end subroutine parmread
-!-----------------------------------------------------------------------------
-! permut.F
-!-----------------------------------------------------------------------------
-      subroutine permut(isym)
-
-      use geometry_data, only: tabperm
-!      use energy_data
-!      use control_data, only:lsecondary
-!      use MD_data
-!      use MPI_data
-!      use map_data
-!      use energy
-!      use geometry
-!      use control
-!      implicit real*8 (a-h,o-z) 
-!      include 'DIMENSIONS'
-!      include 'COMMON.LOCAL'
-!      include 'COMMON.VAR'
-!      include 'COMMON.CHAIN'
-!      include 'COMMON.INTERACT'
-!      include 'COMMON.IOUNITS'
-!      include 'COMMON.GEO'
-!      include 'COMMON.NAMES'
-!      include 'COMMON.CONTROL'
-
-      integer :: n,isym
-!      logical nextp
-!el      external nextp
-      integer,dimension(isym) :: a
-!      parameter(n=symetr)
-!el local variables
-      integer :: kkk,i
-
-      n=isym
-      if (n.eq.1) then
-        tabperm(1,1)=1
-        return
-      endif
-      kkk=0
-      do i=1,n
-      a(i)=i
-      enddo
-   10 print *,(a(i),i=1,n)
-      kkk=kkk+1
-      do i=1,n
-      tabperm(kkk,i)=a(i)
-!      write (iout,*) "tututu", kkk
-      enddo
-      if(nextp(n,a)) go to 10
-      return
-      end subroutine permut
-!-----------------------------------------------------------------------------
-      logical function nextp(n,a)
-
-      integer :: n,i,j,k,t
-!      logical :: nextp
-      integer,dimension(n) :: a
-      i=n-1
-   10 if(a(i).lt.a(i+1)) go to 20
-      i=i-1
-      if(i.eq.0) go to 20
-      go to 10
-   20 j=i+1
-      k=n
-   30 t=a(j)
-      a(j)=a(k)
-      a(k)=t
-      j=j+1
-      k=k-1
-      if(j.lt.k) go to 30
-      j=i
-      if(j.ne.0) go to 40
-      nextp=.false.
-      return
-   40 j=j+1
-      if(a(j).lt.a(i)) go to 40
-      t=a(i)
-      a(i)=a(j)
-      a(j)=t
-      nextp=.true.
-      return
-      end function nextp
+#endif
 !-----------------------------------------------------------------------------
 ! printmat.f
 !-----------------------------------------------------------------------------
@@ -2389,7 +2328,7 @@ write(iout,*) "nloctyp",nloctyp
 !      include 'COMMON.CONTROL'
 !      include 'COMMON.DISTFIT'
 !      include 'COMMON.SETUP'
-      integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity!,&
+      integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift!,ity!,&
 !        ishift_pdb
       logical :: lprn=.true.,fail
       real(kind=8),dimension(3) :: e1,e2,e3
@@ -2421,26 +2360,24 @@ write(iout,*) "nloctyp",nloctyp
       allocate(hfrag(2,maxres/3)) !(2,maxres/3)
       allocate(bfrag(4,maxres/3)) !(4,maxres/3)
 
-!elwrite(iout,*)"poczatek read pdb"
-
       do i=1,100000
         read (ipdbin,'(a80)',end=10) card
 !       write (iout,'(a)') card
         if (card(:5).eq.'HELIX') then
-         nhfrag=nhfrag+1
-         lsecondary=.true.
-         read(card(22:25),*) hfrag(1,nhfrag)
-         read(card(34:37),*) hfrag(2,nhfrag)
+          nhfrag=nhfrag+1
+          lsecondary=.true.
+          read(card(22:25),*) hfrag(1,nhfrag)
+          read(card(34:37),*) hfrag(2,nhfrag)
         endif
         if (card(:5).eq.'SHEET') then
-         nbfrag=nbfrag+1
-         lsecondary=.true.
-         read(card(24:26),*) bfrag(1,nbfrag)
-         read(card(35:37),*) bfrag(2,nbfrag)
+          nbfrag=nbfrag+1
+          lsecondary=.true.
+          read(card(24:26),*) bfrag(1,nbfrag)
+          read(card(35:37),*) bfrag(2,nbfrag)
 !rc----------------------------------------
 !rc  to be corrected !!!
-         bfrag(3,nbfrag)=bfrag(1,nbfrag)
-         bfrag(4,nbfrag)=bfrag(2,nbfrag)
+          bfrag(3,nbfrag)=bfrag(1,nbfrag)
+          bfrag(4,nbfrag)=bfrag(2,nbfrag)
 !rc----------------------------------------
         endif
         if (card(:3).eq.'END') then
@@ -2501,13 +2438,13 @@ write(iout,*) "nloctyp",nloctyp
               endif
               ires=ires-ishift+ishift1
               ires_old=ires
-!              write (iout,*) "ishift",ishift," ires",ires,
-!     &         " ires_old",ires_old
+!              write (iout,*) "ishift",ishift," ires",ires,&
+!               " ires_old",ires_old
               ibeg=0 
             else if (ibeg.eq.2) then
 ! Start a new chain
-!              ishift=-ires_old+ires-1
-!              ishift1=ishift1+1
+              ishift=-ires_old+ires-1 !!!!!
+              ishift1=ishift1-1    !!!!!
 !              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
               ires=ires-ishift+ishift1
               ires_old=ires
@@ -2529,11 +2466,11 @@ write(iout,*) "nloctyp",nloctyp
           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 !            ishift1=ishift1+1
           endif
-!          write (2,*) "ires",ires," res ",res," ity",ity
+!          write (2,*) "ires",ires," res ",res!," ity"!,ity 
           if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
              res.eq.'NHE'.and.atom(:2).eq.'HN') then
             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
-!            write (iout,*) "backbone ",atom 
+!            write (iout,*) "backbone ",atom
 #ifdef DEBUG
             write (iout,'(2i3,2x,a,3f8.3)') &
             ires,itype(ires),res,(c(j,ires),j=1,3)
@@ -2606,32 +2543,32 @@ write(iout,*) "nloctyp",nloctyp
       endif
 !---------------------------------
 !el reallocate tables
-      do i=1,maxres/3
-       do j=1,2
-         hfrag_alloc(j,i)=hfrag(j,i)
-        enddo
-       do j=1,4
-         bfrag_alloc(j,i)=bfrag(j,i)
-        enddo
-      enddo
+!      do i=1,maxres/3
+!      do j=1,2
+!        hfrag_alloc(j,i)=hfrag(j,i)
+!        enddo
+!      do j=1,4
+!        bfrag_alloc(j,i)=bfrag(j,i)
+!        enddo
+!      enddo
 
-      deallocate(hfrag)
-      deallocate(bfrag)
-      allocate(hfrag(2,nres/3)) !(2,maxres/3)
+!      deallocate(hfrag)
+!      deallocate(bfrag)
+!      allocate(hfrag(2,nres/3)) !(2,maxres/3)
 !el      allocate(hfrag(2,nhfrag)) !(2,maxres/3)
 !el      allocate(bfrag(4,nbfrag)) !(4,maxres/3)
-      allocate(bfrag(4,nres/3)) !(4,maxres/3)
+!      allocate(bfrag(4,nres/3)) !(4,maxres/3)
 
-      do i=1,nhfrag
-       do j=1,2
-         hfrag(j,i)=hfrag_alloc(j,i)
-        enddo
-      enddo
-      do i=1,nbfrag
-       do j=1,4
-         bfrag(j,i)=bfrag_alloc(j,i)
-        enddo
-      enddo
+!      do i=1,nhfrag
+!      do j=1,2
+!        hfrag(j,i)=hfrag_alloc(j,i)
+!        enddo
+!      enddo
+!      do i=1,nbfrag
+!      do j=1,4
+!        bfrag(j,i)=bfrag_alloc(j,i)
+!        enddo
+!      enddo
 !el end reallocate tables
 !---------------------------------
       do i=2,nres-1
@@ -2694,13 +2631,52 @@ write(iout,*) "nloctyp",nloctyp
        enddo
       endif
 
-      if(.not.allocated(vbld)) allocate(vbld(2*nres))
-      if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
-      if(.not.allocated(theta)) allocate(theta(nres+2))
+      if(.not.allocated(vbld)) then
+       allocate(vbld(2*nres))
+       do i=1,2*nres
+         vbld(i)=0.d0
+       enddo
+      endif
+      if(.not.allocated(vbld_inv)) then
+       allocate(vbld_inv(2*nres))
+       do i=1,2*nres
+         vbld_inv(i)=0.d0
+       enddo
+      endif
+!!!el
+      if(.not.allocated(theta)) then
+        allocate(theta(nres+2))
+!        allocate(phi(nres+2))
+!        allocate(alph(nres+2))
+!        allocate(omeg(nres+2))
+        do i=1,nres+2
+          theta(i)=0.0d0
+!          phi(i)=0.0d0
+!          alph(i)=0.0d0
+!          omeg(i)=0.0d0
+        enddo
+      endif
+!       allocate(costtab(nres))
+!        allocate(sinttab(nres))
+!        allocate(cost2tab(nres))
+!        allocate(sint2tab(nres))
+!        allocate(xxref(nres))
+!        allocate(yyref(nres))
+!        allocate(zzref(nres)) !(maxres)
+!        do i=1,nres
+!          costtab(i)=0.0d0
+!          sinttab(i)=0.0d0
+!          cost2tab(i)=0.0d0
+!          sint2tab(i)=0.0d0
+!          xxref(i)=0.0d0
+!          yyref(i)=0.0d0
+!          zzref(i)=0.0d0
+!        enddo
+!      endif 
       if(.not.allocated(phi)) allocate(phi(nres+2))
       if(.not.allocated(alph)) allocate(alph(nres+2))
       if(.not.allocated(omeg)) allocate(omeg(nres+2))
-      if(.not.allocated(theta)) allocate(theta(nres+2))
       if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
       if(.not.allocated(phiref)) allocate(phiref(nres+2))
       if(.not.allocated(costtab)) allocate(costtab(nres))
@@ -2710,23 +2686,35 @@ write(iout,*) "nloctyp",nloctyp
       if(.not.allocated(xxref)) allocate(xxref(nres))
       if(.not.allocated(yyref)) allocate(yyref(nres))
       if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
-      if(.not.allocated(theta)) allocate(theta(nres))
-      if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres))
-      if(.not.allocated(theta)) allocate(theta(nres))
+      if(.not.allocated(dc_norm)) then
+!      if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
+        allocate(dc_norm(3,0:2*nres+2))
+        do i=0,2*nres+2
+          dc_norm(1,i)=0.d0
+          dc_norm(2,i)=0.d0
+          dc_norm(3,i)=0.d0
+        enddo
+      endif
       call int_from_cart(.true.,.false.)
       call sc_loc_geom(.true.)
+!      call sc_loc_geom(.false.)
 ! wczesbiej bylo false
       do i=1,nres
         thetaref(i)=theta(i)
         phiref(i)=phi(i)
       enddo
+!      do i=1,2*nres
+!        vbld_inv(i)=0.d0
+!        vbld(i)=0.d0
+!      enddo
       do i=1,nres-1
         do j=1,3
           dc(j,i)=c(j,i+1)-c(j,i)
           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
         enddo
       enddo
-
       do i=2,nres-1
         do j=1,3
           dc(j,i+nres)=c(j,i+nres)-c(j,i)
@@ -2745,9 +2733,9 @@ write(iout,*) "nloctyp",nloctyp
 !        enddo
 !      enddo
 !
-      allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
-      allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
-      allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
+      if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
+      if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
+      if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
 !-----------------------------
       kkk=1
       lll=0
@@ -2850,29 +2838,10 @@ write(iout,*) "nloctyp",nloctyp
         enddo
       enddo
       ishift_pdb=ishift
-!---------------------
-! el reallocate array
-      do i=1,2*nres+2
-        do kkk=1,nperm
-          cref_alloc(1,i,kkk)=cref(1,i,kkk)
-          cref_alloc(2,i,kkk)=cref(2,i,kkk)
-          cref_alloc(3,i,kkk)=cref(3,i,kkk)
-        enddo
-      enddo
-!el      deallocate(cref)
-!el      allocate(cref(3,2*nres+2,nperm)) !(3,maxres2+2,maxperm)
-
-      do i=1,2*nres+2
-        do kkk=1,nperm
-          cref(1,i,kkk)=cref_alloc(1,i,kkk)
-          cref(2,i,kkk)=cref_alloc(2,i,kkk)
-          cref(3,i,kkk)=cref_alloc(3,i,kkk)
-        enddo
-      enddo
-!---------------------
+
       return
       end subroutine readpdb
-#ifndef WHAM_RUN
+#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
 !-----------------------------------------------------------------------------
 ! readrtns_CSA.F
 !-----------------------------------------------------------------------------
@@ -3603,7 +3572,6 @@ write(iout,*) "nloctyp",nloctyp
       use energy_data, only: usampl
       use csa_data
       use MPI_data
-!      use MD
       use control_data, only:out1file
       use control, only: getenv_loc
 !      implicit real*8 (a-h,o-z)
@@ -3939,8 +3907,12 @@ write(iout,*) "nloctyp",nloctyp
         thetname(:ilen(thetname))
       write (iout,*) "Rotamer parameter file          : ",&
         rotname(:ilen(rotname))
+!el----
+#ifndef CRYST_THETA
       write (iout,*) "Thetpdb parameter file          : ",&
         thetname_pdb(:ilen(thetname_pdb))
+#endif
+!el
       write (iout,*) "Threading database              : ",&
         patname(:ilen(patname))
       if (lentmp.ne.0) &