Merge branch 'UCGM' of mmka.chem.univ.gda.pl:unres4 into UCGM

[unres4.git] / source / unres / io_config.f90

diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90
index 91daac1..87ddf94 100644 (file)
--- a/source/unres/io_config.f90
+++ b/source/unres/io_config.f90
@@ -478,7 +478,7 @@
 
 !     write (iout,100)
 !      do i=1,nres
 
 !     write (iout,100)
 !      do i=1,nres

-!        write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
+!        write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&

 !          c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
 !      enddo
 !  100 format (//'              alpha-carbon coordinates       ',&
 !          c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
 !      enddo
 !  100 format (//'              alpha-carbon coordinates       ',&

@@ -529,9 +529,9 @@
 ! read secondary structure prediction from JPRED here!
 !      read(isecpred,'(A80)',err=100,end=100) line
 !      read(line,'(f10.3)',err=110) ftors
 ! read secondary structure prediction from JPRED here!
 !      read(isecpred,'(A80)',err=100,end=100) line
 !      read(line,'(f10.3)',err=110) ftors

-       read(isecpred,'(f10.3)',err=110) ftors
+       read(isecpred,'(f10.3)',err=110) ftors(1)

 
 

-      write (iout,*) 'FTORS factor =',ftors
+      write (iout,*) 'FTORS factor =',ftors(1)

 ! initialize secstruc to any
        do i=1,nres
         secstruc(i) ='-'
 ! initialize secstruc to any
        do i=1,nres
         secstruc(i) ='-'

@@ -553,6 +553,7 @@
       
       ii=0
       do i=1,nres
       
       ii=0
       do i=1,nres

+         ftors(i)=ftors(1)

         if ( secstruc(i) .eq. 'H') then
 ! Helix restraints for this residue               
            ii=ii+1 
         if ( secstruc(i) .eq. 'H') then
 ! Helix restraints for this residue               
            ii=ii+1 

@@ -735,11 +736,11 @@
       character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
       logical :: lprint,LaTeX
       real(kind=8),dimension(3,3,maxlob) :: blower     !(3,3,maxlob)
       character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
       logical :: lprint,LaTeX
       real(kind=8),dimension(3,3,maxlob) :: blower     !(3,3,maxlob)

-      real(kind=8),dimension(13) :: b
+      real(kind=8),dimension(13) :: buse

       character(len=3) :: lancuch      !,ucase
 !el  local variables
       integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
       character(len=3) :: lancuch      !,ucase
 !el  local variables
       integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm

-      integer :: maxinter,junk,kk,ii
+      integer :: maxinter,junk,kk,ii,ncatprotparm

       real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
                 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
                 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
       real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
                 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
                 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&

@@ -776,12 +777,12 @@
 !
       allocate(dsc(ntyp1)) !(ntyp1)
       allocate(dsc_inv(ntyp1)) !(ntyp1)
 !
       allocate(dsc(ntyp1)) !(ntyp1)
       allocate(dsc_inv(ntyp1)) !(ntyp1)

+      allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
+      allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
+      allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)

       allocate(nbondterm(ntyp)) !(ntyp)
       allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
       allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
       allocate(nbondterm(ntyp)) !(ntyp)
       allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
       allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)

-      allocate(msc(ntyp+1)) !(ntyp+1)
-      allocate(isc(ntyp+1)) !(ntyp+1)
-      allocate(restok(ntyp+1)) !(ntyp+1)

       allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
       allocate(long_r_sidechain(ntyp))
       allocate(short_r_sidechain(ntyp))
       allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
       allocate(long_r_sidechain(ntyp))
       allocate(short_r_sidechain(ntyp))

@@ -789,6 +790,10 @@
       dsc_inv(:)=0.0d0
 
 #ifdef CRYST_BOND
       dsc_inv(:)=0.0d0
 
 #ifdef CRYST_BOND

+      allocate(msc(ntyp+1)) !(ntyp+1)
+      allocate(isc(ntyp+1)) !(ntyp+1)
+      allocate(restok(ntyp+1)) !(ntyp+1)
+

       read (ibond,*) vbldp0,akp,mp,ip,pstok
       do i=1,ntyp
         nbondterm(i)=1
       read (ibond,*) vbldp0,akp,mp,ip,pstok
       do i=1,ntyp
         nbondterm(i)=1

@@ -801,10 +806,19 @@
         endif
       enddo
 #else
         endif
       enddo
 #else

-      read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp,ip,pstok
-      do i=1,ntyp
+      mp(:)=0.0d0
+      ip(:)=0.0d0
+      msc(:,:)=0.0d0
+      isc(:,:)=0.0d0
+
+      allocate(msc(ntyp+1,5)) !(ntyp+1)
+      allocate(isc(ntyp+1,5)) !(ntyp+1)
+      allocate(restok(ntyp+1,5)) !(ntyp+1)
+
+      read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
+      do i=1,ntyp_molec(1)

         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&

-         j=1,nbondterm(i)),msc(i),isc(i),restok(i)
+         j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)

         dsc(i) = vbldsc0(1,i)
         if (i.eq.10) then
           dsc_inv(i)=0.0D0
         dsc(i) = vbldsc0(1,i)
         if (i.eq.10) then
           dsc_inv(i)=0.0D0

@@ -813,20 +827,57 @@
         endif
       enddo
 #endif
         endif
       enddo
 #endif

+      read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
+      do i=1,ntyp_molec(2)
+        nbondterm_nucl(i)=1
+        read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
+!        dsc(i) = vbldsc0_nucl(1,i)
+!        if (i.eq.10) then
+!          dsc_inv(i)=0.0D0
+!        else
+!          dsc_inv(i)=1.0D0/dsc(i)
+!        endif
+      enddo
+!      read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
+!      do i=1,ntyp_molec(2)
+!        read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),& 
+!        aksc_nucl(j,i),abond0_nucl(j,i),&
+!         j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
+!        dsc(i) = vbldsc0(1,i)
+!        if (i.eq.10) then
+!          dsc_inv(i)=0.0D0
+!        else
+!          dsc_inv(i)=1.0D0/dsc(i)
+!        endif
+!      enddo
+

       if (lprint) then
         write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
          'inertia','Pstok'
       if (lprint) then
         write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
          'inertia','Pstok'

-        write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
+        write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)

         do i=1,ntyp
         do i=1,ntyp

-          write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
-            vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
+          write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
+            vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)

           do j=2,nbondterm(i)
             write (iout,'(13x,3f10.5)') &
               vbldsc0(j,i),aksc(j,i),abond0(j,i)
           enddo
         enddo
       endif
           do j=2,nbondterm(i)
             write (iout,'(13x,3f10.5)') &
               vbldsc0(j,i),aksc(j,i),abond0(j,i)
           enddo
         enddo
       endif

+            do i=1,ntyp_molec(5)
+             read(iion,*) msc(i,5),restok(i,5)
+             print *,msc(i,5),restok(i,5)
+            enddo
+            ip(5)=0.2
+!            isc(5)=0.2
+            read (iion,*) ncatprotparm
+            allocate(catprm(ncatprotparm,4))
+            do k=1,4
+            read (iion,*)  (catprm(i,k),i=1,ncatprotparm)
+            enddo
+            print *, catprm
+!            read (iion,*) (vcatprm(k),k=1,ncatprotpram)

 !----------------------------------------------------
       allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
       allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp))     !(-ntyp:ntyp)
 !----------------------------------------------------
       allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
       allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp))     !(-ntyp:ntyp)

@@ -912,7 +963,7 @@
        '    ATHETA0   ','         A1   ','        A2    ',&
        '        B1    ','         B2   '        
         do i=1,ntyp
        '    ATHETA0   ','         A1   ','        A2    ',&
        '        B1    ','         B2   '        
         do i=1,ntyp

-          write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
+          write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&

               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-))') &
               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-))') &

@@ -920,7 +971,7 @@
        '     ALPH0    ','      ALPH1   ','     ALPH2    ',&
        '     ALPH3    ','    SIGMA0C   '        
         do i=1,ntyp
        '     ALPH0    ','      ALPH1   ','     ALPH2    ',&
        '     ALPH3    ','    SIGMA0C   '        
         do i=1,ntyp

-          write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
+          write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&

             (polthet(j,i),j=0,3),sigc0(i) 
         enddo
         write (iout,'(/a/9x,5a/79(1h-))') &
             (polthet(j,i),j=0,3),sigc0(i) 
         enddo
         write (iout,'(/a/9x,5a/79(1h-))') &

@@ -928,7 +979,7 @@
        '    THETA0    ','     SIGMA0   ','        G1    ',&
        '        G2    ','         G3   '        
         do i=1,ntyp
        '    THETA0    ','     SIGMA0   ','        G1    ',&
        '        G2    ','         G3   '        
         do i=1,ntyp

-          write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),&
+          write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&

              sig0(i),(gthet(j,i),j=1,3)
         enddo
        else
              sig0(i),(gthet(j,i),j=1,3)
         enddo
        else

@@ -939,7 +990,7 @@
        '     theta0   ','    a1*10^2   ','   a2*10^2    ',&
        '   b1*10^1    ','    b2*10^1   '        
         do i=1,ntyp
        '     theta0   ','    a1*10^2   ','   a2*10^2    ',&
        '   b1*10^1    ','    b2*10^1   '        
         do i=1,ntyp

-          write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
+          write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&

               a0thet(i),(100*athet(j,i,1,1),j=1,2),&
               (10*bthet(j,i,1,1),j=1,2)
         enddo
               a0thet(i),(100*athet(j,i,1,1),j=1,2),&
               (10*bthet(j,i,1,1),j=1,2)
         enddo

@@ -948,7 +999,7 @@
        ' alpha0       ','  alph1       ',' alph2        ',&
        ' alhp3        ','   sigma0c    '        
        do i=1,ntyp
        ' alpha0       ','  alph1       ',' alph2        ',&
        ' alhp3        ','   sigma0c    '        
        do i=1,ntyp

-          write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
+          write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&

             (polthet(j,i),j=0,3),sigc0(i) 
        enddo
        write (iout,'(/a/9x,5a/79(1h-))') &
             (polthet(j,i),j=0,3),sigc0(i) 
        enddo
        write (iout,'(/a/9x,5a/79(1h-))') &

@@ -956,7 +1007,7 @@
        '    theta0    ','  sigma0*10^2 ','      G1*10^-1',&
        '        G2    ','   G3*10^1    '        
        do i=1,ntyp
        '    theta0    ','  sigma0*10^2 ','      G1*10^-1',&
        '        G2    ','   G3*10^1    '        
        do i=1,ntyp

-          write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
+          write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&

              100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
        enddo
       endif
              100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
        enddo
       endif

@@ -1207,6 +1258,65 @@
       close (ithep_pdb)
 #endif
       close(ithep)
       close (ithep_pdb)
 #endif
       close(ithep)

+!--------------- Reading theta parameters for nucleic acid-------
+      read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
+      ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
+      nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
+      allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
+      allocate(aa0thet_nucl(nthetyp_nucl+1,&
+        nthetyp_nucl+1,nthetyp_nucl+1))
+!(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+      allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
+        nthetyp_nucl+1,nthetyp_nucl+1))
+!(maxtheterm,-maxthetyp1:maxthetyp1,&
+!        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+      allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+        nthetyp_nucl+1,nthetyp_nucl+1))
+      allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+        nthetyp_nucl+1,nthetyp_nucl+1))
+      allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+        nthetyp_nucl+1,nthetyp_nucl+1))
+      allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+        nthetyp_nucl+1,nthetyp_nucl+1))
+!(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
+!        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+      allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
+         nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
+      allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
+         nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
+
+!(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
+!        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+
+      read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
+
+      aa0thet_nucl(:,:,:)=0.0d0
+      aathet_nucl(:,:,:,:)=0.0d0
+      bbthet_nucl(:,:,:,:,:)=0.0d0
+      ccthet_nucl(:,:,:,:,:)=0.0d0
+      ddthet_nucl(:,:,:,:,:)=0.0d0
+      eethet_nucl(:,:,:,:,:)=0.0d0
+      ffthet_nucl(:,:,:,:,:,:)=0.0d0
+      ggthet_nucl(:,:,:,:,:,:)=0.0d0
+
+      do i=1,nthetyp_nucl
+        do j=1,nthetyp_nucl
+          do k=1,nthetyp_nucl
+            read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
+            read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
+            read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
+            read (ithep_nucl,*,end=111,err=111) &
+            (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
+            (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
+            (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
+            (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
+            read (ithep_nucl,*,end=111,err=111) &
+           (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
+              ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
+              llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
+          enddo
+        enddo
+      enddo

 
 !-------------------------------------------
       allocate(nlob(ntyp1)) !(ntyp1)
 
 !-------------------------------------------
       allocate(nlob(ntyp1)) !(ntyp1)

@@ -1285,7 +1395,7 @@
          nlobi=nlob(i)
           if (nlobi.gt.0) then
             if (LaTeX) then
          nlobi=nlob(i)
           if (nlobi.gt.0) then
             if (LaTeX) then

-              write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
+              write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&

                ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
                                    'log h',(bsc(j,i),j=1,nlobi)
                ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
                                    'log h',(bsc(j,i),j=1,nlobi)

@@ -1326,6 +1436,24 @@
          enddo  
        endif
       enddo
          enddo  
        endif
       enddo

+!---------reading nucleic acid parameters for rotamers-------------------
+      allocate(sc_parmin_nucl(9,ntyp_molec(2)))      !(maxsccoef,ntyp)
+      do i=1,ntyp_molec(2)
+        read (irotam_nucl,*,end=112,err=112)
+        do j=1,9
+          read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
+        enddo
+      enddo
+      close(irotam_nucl)
+      if (lprint) then
+        write (iout,*)
+        write (iout,*) "Base rotamer parameters"
+        do i=1,ntyp_molec(2)
+          write (iout,'(a)') restyp(i,2)
+          write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
+        enddo
+      endif
+

 !
 ! Read the parameters of the probability distribution/energy expression
 ! of the side chains.
 !
 ! Read the parameters of the probability distribution/energy expression
 ! of the side chains.

@@ -1595,6 +1723,48 @@
       enddo
       endif
 #endif
       enddo
       endif
 #endif

+      allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
+      read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
+!      print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
+!el from energy module---------
+      allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
+      allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
+
+      allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
+      allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
+      allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
+      allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
+
+      allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
+      allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
+!el---------------------------
+      nterm_nucl(:,:)=0
+      nlor_nucl(:,:)=0
+!el--------------------
+      read (itorp_nucl,*,end=113,err=113) &
+        (itortyp_nucl(i),i=1,ntyp_molec(2))
+!        print *,itortyp_nucl(:)
+!c      write (iout,*) 'ntortyp',ntortyp
+      do i=1,ntortyp_nucl
+        do j=1,ntortyp_nucl
+          read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
+!           print *,nterm_nucl(i,j),nlor_nucl(i,j)
+          v0ij=0.0d0
+          si=-1.0d0
+          do k=1,nterm_nucl(i,j)
+            read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
+            v0ij=v0ij+si*v1_nucl(k,i,j)
+            si=-si
+          enddo
+          do k=1,nlor_nucl(i,j)
+            read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
+              vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
+            v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
+          enddo
+          v0_nucl(i,j)=v0ij
+        enddo
+      enddo
+

 ! Read of Side-chain backbone correlation parameters
 ! Modified 11 May 2012 by Adasko
 !CC
 ! Read of Side-chain backbone correlation parameters
 ! Modified 11 May 2012 by Adasko
 !CC

@@ -1738,6 +1908,7 @@
 
 #endif      
       if (lprint) then
 
 #endif      
       if (lprint) then

+        l=3

         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

@@ -1782,88 +1953,88 @@
 
       do i=0,nloctyp-1
         read (ifourier,*,end=115,err=115)
 
       do i=0,nloctyp-1
         read (ifourier,*,end=115,err=115)

-        read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
+        read (ifourier,*,end=115,err=115) (buse(ii),ii=1,13)

         if (lprint) then
           write (iout,*) 'Type',i
         if (lprint) then
           write (iout,*) 'Type',i

-          write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
+          write (iout,'(a,i2,a,f10.5)') ('buse(',ii,')=',buse(ii),ii=1,13)

         endif
         endif

-        B1(1,i)  = b(3)
-        B1(2,i)  = b(5)
-        B1(1,-i) = b(3)
-        B1(2,-i) = -b(5)
-!        b1(1,i)=0.0d0
-!        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(1,i)=0.0d0
-!        b1tilde(2,i)=0.0d0
-        B2(1,i)  = b(2)
-        B2(2,i)  = b(4)
-        B2(1,-i)  =b(2)
-        B2(2,-i)  =-b(4)
-
-!        b2(1,i)=0.0d0
-!        b2(2,i)=0.0d0
-        CC(1,1,i)= b(7)
-        CC(2,2,i)=-b(7)
-        CC(2,1,i)= b(9)
-        CC(1,2,i)= b(9)
-        CC(1,1,-i)= b(7)
-        CC(2,2,-i)=-b(7)
-        CC(2,1,-i)=-b(9)
-        CC(1,2,-i)=-b(9)
+        B1(1,i)  = buse(3)
+        B1(2,i)  = buse(5)
+        B1(1,-i) = buse(3)
+        B1(2,-i) = -buse(5)
+!        buse1(1,i)=0.0d0
+!        buse1(2,i)=0.0d0
+        B1tilde(1,i) = buse(3)
+        B1tilde(2,i) =-buse(5)
+        B1tilde(1,-i) =-buse(3)
+        B1tilde(2,-i) =buse(5)
+!        buse1tilde(1,i)=0.0d0
+!        buse1tilde(2,i)=0.0d0
+        B2(1,i)  = buse(2)
+        B2(2,i)  = buse(4)
+        B2(1,-i)  =buse(2)
+        B2(2,-i)  =-buse(4)
+
+!        buse2(1,i)=0.0d0
+!        buse2(2,i)=0.0d0
+        CC(1,1,i)= buse(7)
+        CC(2,2,i)=-buse(7)
+        CC(2,1,i)= buse(9)
+        CC(1,2,i)= buse(9)
+        CC(1,1,-i)= buse(7)
+        CC(2,2,-i)=-buse(7)
+        CC(2,1,-i)=-buse(9)
+        CC(1,2,-i)=-buse(9)

 !        CC(1,1,i)=0.0d0
 !        CC(2,2,i)=0.0d0
 !        CC(2,1,i)=0.0d0
 !        CC(1,2,i)=0.0d0
 !        CC(1,1,i)=0.0d0
 !        CC(2,2,i)=0.0d0
 !        CC(2,1,i)=0.0d0
 !        CC(1,2,i)=0.0d0

-        Ctilde(1,1,i)=b(7)
-        Ctilde(1,2,i)=b(9)
-        Ctilde(2,1,i)=-b(9)
-        Ctilde(2,2,i)=b(7)
-        Ctilde(1,1,-i)=b(7)
-        Ctilde(1,2,-i)=-b(9)
-        Ctilde(2,1,-i)=b(9)
-        Ctilde(2,2,-i)=b(7)
+        Ctilde(1,1,i)=buse(7)
+        Ctilde(1,2,i)=buse(9)
+        Ctilde(2,1,i)=-buse(9)
+        Ctilde(2,2,i)=buse(7)
+        Ctilde(1,1,-i)=buse(7)
+        Ctilde(1,2,-i)=-buse(9)
+        Ctilde(2,1,-i)=buse(9)
+        Ctilde(2,2,-i)=buse(7)

 
 !        Ctilde(1,1,i)=0.0d0
 !        Ctilde(1,2,i)=0.0d0
 !        Ctilde(2,1,i)=0.0d0
 !        Ctilde(2,2,i)=0.0d0
 
 !        Ctilde(1,1,i)=0.0d0
 !        Ctilde(1,2,i)=0.0d0
 !        Ctilde(2,1,i)=0.0d0
 !        Ctilde(2,2,i)=0.0d0

-        DD(1,1,i)= b(6)
-        DD(2,2,i)=-b(6)
-        DD(2,1,i)= b(8)
-        DD(1,2,i)= b(8)
-        DD(1,1,-i)= b(6)
-        DD(2,2,-i)=-b(6)
-        DD(2,1,-i)=-b(8)
-        DD(1,2,-i)=-b(8)
+        DD(1,1,i)= buse(6)
+        DD(2,2,i)=-buse(6)
+        DD(2,1,i)= buse(8)
+        DD(1,2,i)= buse(8)
+        DD(1,1,-i)= buse(6)
+        DD(2,2,-i)=-buse(6)
+        DD(2,1,-i)=-buse(8)
+        DD(1,2,-i)=-buse(8)

 !        DD(1,1,i)=0.0d0
 !        DD(2,2,i)=0.0d0
 !        DD(2,1,i)=0.0d0
 !        DD(1,2,i)=0.0d0
 !        DD(1,1,i)=0.0d0
 !        DD(2,2,i)=0.0d0
 !        DD(2,1,i)=0.0d0
 !        DD(1,2,i)=0.0d0

-        Dtilde(1,1,i)=b(6)
-        Dtilde(1,2,i)=b(8)
-        Dtilde(2,1,i)=-b(8)
-        Dtilde(2,2,i)=b(6)
-        Dtilde(1,1,-i)=b(6)
-        Dtilde(1,2,-i)=-b(8)
-        Dtilde(2,1,-i)=b(8)
-        Dtilde(2,2,-i)=b(6)
+        Dtilde(1,1,i)=buse(6)
+        Dtilde(1,2,i)=buse(8)
+        Dtilde(2,1,i)=-buse(8)
+        Dtilde(2,2,i)=buse(6)
+        Dtilde(1,1,-i)=buse(6)
+        Dtilde(1,2,-i)=-buse(8)
+        Dtilde(2,1,-i)=buse(8)
+        Dtilde(2,2,-i)=buse(6)

 
 !        Dtilde(1,1,i)=0.0d0
 !        Dtilde(1,2,i)=0.0d0
 !        Dtilde(2,1,i)=0.0d0
 !        Dtilde(2,2,i)=0.0d0
 
 !        Dtilde(1,1,i)=0.0d0
 !        Dtilde(1,2,i)=0.0d0
 !        Dtilde(2,1,i)=0.0d0
 !        Dtilde(2,2,i)=0.0d0

-        EE(1,1,i)= b(10)+b(11)
-        EE(2,2,i)=-b(10)+b(11)
-        EE(2,1,i)= b(12)-b(13)
-        EE(1,2,i)= b(12)+b(13)
-        EE(1,1,-i)= b(10)+b(11)
-        EE(2,2,-i)=-b(10)+b(11)
-        EE(2,1,-i)=-b(12)+b(13)
-        EE(1,2,-i)=-b(12)-b(13)
+        EE(1,1,i)= buse(10)+buse(11)
+        EE(2,2,i)=-buse(10)+buse(11)
+        EE(2,1,i)= buse(12)-buse(13)
+        EE(1,2,i)= buse(12)+buse(13)
+        EE(1,1,-i)= buse(10)+buse(11)
+        EE(2,2,-i)=-buse(10)+buse(11)
+        EE(2,1,-i)=-buse(12)+buse(13)
+        EE(1,2,-i)=-buse(12)-buse(13)

 
 !        ee(1,1,i)=1.0d0
 !        ee(2,2,i)=1.0d0
 
 !        ee(1,1,i)=1.0d0
 !        ee(2,2,i)=1.0d0

@@ -1927,6 +2098,7 @@
       allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
       allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
       allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
       allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
       allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
       allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)

+

       allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
       allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
       allocate(epslip(ntyp,ntyp))
       allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
       allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
       allocate(epslip(ntyp,ntyp))

@@ -1965,7 +2137,7 @@
          call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
          write (iout,'(/a)') 'One-body parameters:'
          write (iout,'(a,4x,a)') 'residue','sigma'
          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)
+         write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)

         endif
 !      goto 50
 !----------------------- LJK potential --------------------------------
         endif
 !      goto 50
 !----------------------- LJK potential --------------------------------

@@ -1979,13 +2151,15 @@
          call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
          write (iout,'(/a)') 'One-body parameters:'
          write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
          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),&
+          write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&

                 i=1,ntyp)
         endif
 !      goto 50
 !---------------------- GB or BP potential -----------------------------
        case(3:4)
 !   30 do i=1,ntyp
                 i=1,ntyp)
         endif
 !      goto 50
 !---------------------- GB or BP potential -----------------------------
        case(3:4)
 !   30 do i=1,ntyp

+!        print *,"I AM in SCELE",scelemode
+        if (scelemode.eq.0) then

         do i=1,ntyp
          read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
         enddo
         do i=1,ntyp
          read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
         enddo

@@ -2010,9 +2184,84 @@
          write (iout,'(/a)') 'One-body parameters:'
          write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',&
                '    chip  ','    alph  '
          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),&
+         write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&

                              chip(i),alp(i),i=1,ntyp)
         endif
                              chip(i),alp(i),i=1,ntyp)
         endif

+       else
+!      print *,ntyp,"NTYP"
+      allocate(icharge(ntyp1))
+!      print *,ntyp,icharge(i)
+      icharge(:)=0
+      read (isidep,*) (icharge(i),i=1,ntyp)
+      print *,ntyp,icharge(i)
+!      if(.not.allocated(eps)) allocate(eps(-ntyp
+!c      write (2,*) "icharge",(icharge(i),i=1,ntyp)
+       allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
+       allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
+       allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
+       allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
+       allocate(epsintab(ntyp,ntyp))
+       allocate(dtail(2,ntyp,ntyp))
+       allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
+       allocate(wqdip(2,ntyp,ntyp))
+       allocate(wstate(4,ntyp,ntyp))
+       allocate(dhead(2,2,ntyp,ntyp))
+       allocate(nstate(ntyp,ntyp))
+      if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
+      if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
+      do i=1,ntyp
+       do j=1,i
+!        write (*,*) "Im in ALAB", i, " ", j
+        read(isidep,*) &
+       eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
+       (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
+       chis(i,j),chis(j,i), &
+       nstate(i,j),(wstate(k,i,j),k=1,4), &
+       dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
+       dtail(1,i,j),dtail(2,i,j), &
+       epshead(i,j),sig0head(i,j), &
+       rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
+       alphapol(i,j),alphapol(j,i), &
+       (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
+!       print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i) 
+       END DO
+      END DO
+      DO i = 1, ntyp
+       DO j = i+1, ntyp
+        eps(i,j) = eps(j,i)
+        sigma(i,j) = sigma(j,i)
+        sigmap1(i,j)=sigmap1(j,i)
+        sigmap2(i,j)=sigmap2(j,i)
+        sigiso1(i,j)=sigiso1(j,i)
+        sigiso2(i,j)=sigiso2(j,i)
+!        print *,"ATU",sigma(j,i),sigma(i,j),i,j
+        nstate(i,j) = nstate(j,i)
+        dtail(1,i,j) = dtail(1,j,i)
+        dtail(2,i,j) = dtail(2,j,i)
+        DO k = 1, 4
+         alphasur(k,i,j) = alphasur(k,j,i)
+         wstate(k,i,j) = wstate(k,j,i)
+         alphiso(k,i,j) = alphiso(k,j,i)
+        END DO
+
+        dhead(2,1,i,j) = dhead(1,1,j,i)
+        dhead(2,2,i,j) = dhead(1,2,j,i)
+        dhead(1,1,i,j) = dhead(2,1,j,i)
+        dhead(1,2,i,j) = dhead(2,2,j,i)
+
+        epshead(i,j) = epshead(j,i)
+        sig0head(i,j) = sig0head(j,i)
+
+        DO k = 1, 2
+         wqdip(k,i,j) = wqdip(k,j,i)
+        END DO
+
+        wquad(i,j) = wquad(j,i)
+        epsintab(i,j) = epsintab(j,i)
+!        if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
+       END DO
+      END DO
+      endif

 !      goto 50
 !--------------------- GBV potential -----------------------------------
        case(5)
 !      goto 50
 !--------------------- GBV potential -----------------------------------
        case(5)

@@ -2027,7 +2276,7 @@
           write (iout,'(/a)') 'One-body parameters:'
           write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',&
               's||/s_|_^2','    chip  ','    alph  '
           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),&
+          write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&

                    sigii(i),chip(i),alp(i),i=1,ntyp)
         endif
        case default
                    sigii(i),chip(i),alp(i),i=1,ntyp)
         endif
        case default

@@ -2036,13 +2285,16 @@
       end select
       continue
       close (isidep)
       end select
       continue
       close (isidep)

+

 !-----------------------------------------------------------------------
 ! Calculate the "working" parameters of SC interactions.
 
 !el from module energy - COMMON.INTERACT-------
 !-----------------------------------------------------------------------
 ! Calculate the "working" parameters of SC interactions.
 
 !el from module energy - COMMON.INTERACT-------

-      allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
+      allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
+      if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)

       allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
       allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)

-      allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
+      if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
+      allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)

       allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
         dcavtub(ntyp1))
       allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
       allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
         dcavtub(ntyp1))
       allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&

@@ -2051,9 +2303,11 @@
       bb_aq(:,:)=0.0D0
       aa_lip(:,:)=0.0D0
       bb_lip(:,:)=0.0D0
       bb_aq(:,:)=0.0D0
       aa_lip(:,:)=0.0D0
       bb_lip(:,:)=0.0D0

+         if (scelemode.eq.0) then

       chi(:,:)=0.0D0
       sigma(:,:)=0.0D0
       r0(:,:)=0.0D0
       chi(:,:)=0.0D0
       sigma(:,:)=0.0D0
       r0(:,:)=0.0D0

+        endif

       acavtub(:)=0.0d0
       bcavtub(:)=0.0d0
       ccavtub(:)=0.0d0
       acavtub(:)=0.0d0
       bcavtub(:)=0.0d0
       ccavtub(:)=0.0d0

@@ -2069,6 +2323,7 @@
           epslip(i,j)=epslip(j,i)
         enddo
       enddo
           epslip(i,j)=epslip(j,i)
         enddo
       enddo

+         if (scelemode.eq.0) then

       do i=1,ntyp
         do j=i,ntyp
           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
       do i=1,ntyp
         do j=i,ntyp
           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)

@@ -2077,6 +2332,7 @@
           rs0(j,i)=rs0(i,j)
         enddo
       enddo
           rs0(j,i)=rs0(i,j)
         enddo
       enddo

+      endif

       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
        'Working parameters of the SC interactions:',&
        '     a    ','     b    ','   augm   ','  sigma ','   r0   ',&
       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
        'Working parameters of the SC interactions:',&
        '     a    ','     b    ','   augm   ','  sigma ','   r0   ',&

@@ -2086,6 +2342,7 @@
          epsij=eps(i,j)
          if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
            rrij=sigma(i,j)
          epsij=eps(i,j)
          if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
            rrij=sigma(i,j)

+!            print *,"SIGMA ZLA?",sigma(i,j)

           else
            rrij=rr0(i)+rr0(j)
           endif
           else
            rrij=rr0(i)+rr0(j)
           endif

@@ -2096,6 +2353,7 @@
          sigeps=dsign(1.0D0,epsij)
          epsij=dabs(epsij)
          aa_aq(i,j)=epsij*rrij*rrij
          sigeps=dsign(1.0D0,epsij)
          epsij=dabs(epsij)
          aa_aq(i,j)=epsij*rrij*rrij

+          print *,"ADASKO",epsij,rrij,expon

          bb_aq(i,j)=-sigeps*epsij*rrij
          aa_aq(j,i)=aa_aq(i,j)
          bb_aq(j,i)=bb_aq(i,j)
          bb_aq(i,j)=-sigeps*epsij*rrij
          aa_aq(j,i)=aa_aq(i,j)
          bb_aq(j,i)=bb_aq(i,j)

@@ -2107,7 +2365,7 @@
           aa_lip(j,i)=aa_lip(i,j)
           bb_lip(j,i)=bb_lip(i,j)
 !C          write(iout,*) aa_lip
           aa_lip(j,i)=aa_lip(i,j)
           bb_lip(j,i)=bb_lip(i,j)
 !C          write(iout,*) aa_lip

-         if (ipot.gt.2) then
+         if ((ipot.gt.2).and. (scelemode.eq.0)) then

            sigt1sq=sigma0(i)**2
            sigt2sq=sigma0(j)**2
            sigii1=sigii(i)
            sigt1sq=sigma0(i)**2
            sigt2sq=sigma0(j)**2
            sigii1=sigii(i)

@@ -2139,11 +2397,90 @@
           endif
          if (lprint) then
             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
           endif
          if (lprint) then
             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &

-            restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
+            restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&

             sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
          endif
         enddo
       enddo
             sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
          endif
         enddo
       enddo

+
+      allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
+      allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
+      allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
+      allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+      allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
+      allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+      allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
+      allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+      allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+      allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+      allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+      allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+      allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+      allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+      allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
+      allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
+
+!      augm(:,:)=0.0D0
+!      chip(:)=0.0D0
+!      alp(:)=0.0D0
+!      sigma0(:)=0.0D0
+!      sigii(:)=0.0D0
+!      rr0(:)=0.0D0
+   
+      read (isidep_nucl,*) ipot_nucl
+!      print *,"TU?!",ipot_nucl
+      if (ipot_nucl.eq.1) then
+        do i=1,ntyp_molec(2)
+          do j=i,ntyp_molec(2)
+            read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
+            elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
+          enddo
+        enddo
+      else
+        do i=1,ntyp_molec(2)
+          do j=i,ntyp_molec(2)
+            read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
+               chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
+               elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
+          enddo
+        enddo
+      endif
+!      rpp(1,1)=2**(1.0/6.0)*5.16158
+      do i=1,ntyp_molec(2)
+        do j=i,ntyp_molec(2)
+          rrij=sigma_nucl(i,j)
+          r0_nucl(i,j)=rrij
+          r0_nucl(j,i)=rrij
+          rrij=rrij**expon
+          epsij=4*eps_nucl(i,j)
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_nucl(i,j)=epsij*rrij*rrij
+          bb_nucl(i,j)=-sigeps*epsij*rrij
+          ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
+          ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
+          ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
+          ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
+          sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
+         2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
+        enddo
+        do j=1,i-1
+          aa_nucl(i,j)=aa_nucl(j,i)
+          bb_nucl(i,j)=bb_nucl(j,i)
+          ael3_nucl(i,j)=ael3_nucl(j,i)
+          ael6_nucl(i,j)=ael6_nucl(j,i)
+          ael63_nucl(i,j)=ael63_nucl(j,i)
+          ael32_nucl(i,j)=ael32_nucl(j,i)
+          elpp3_nucl(i,j)=elpp3_nucl(j,i)
+          elpp6_nucl(i,j)=elpp6_nucl(j,i)
+          elpp63_nucl(i,j)=elpp63_nucl(j,i)
+          elpp32_nucl(i,j)=elpp32_nucl(j,i)
+          eps_nucl(i,j)=eps_nucl(j,i)
+          sigma_nucl(i,j)=sigma_nucl(j,i)
+          sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
+        enddo
+      enddo
+

       write(iout,*) "tube param"
       read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
       ccavtubpep,dcavtubpep,tubetranenepep
       write(iout,*) "tube param"
       read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
       ccavtubpep,dcavtubpep,tubetranenepep

@@ -2163,6 +2500,174 @@
        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
        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

+!-----------------READING SC BASE POTENTIALS-----------------------------
+      allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))      
+      allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
+      allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
+      allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
+      allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
+      allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
+      allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
+      allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
+      allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
+      allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
+      allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
+      allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
+      allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
+      allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
+      allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
+      allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
+
+
+      do i=1,ntyp_molec(1)
+       do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
+        write (*,*) "Im in ", i, " ", j
+        read(isidep_scbase,*) &
+        eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
+        chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
+         write(*,*) "eps",eps_scbase(i,j)
+        read(isidep_scbase,*) &
+       (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
+       chis_scbase(i,j,1),chis_scbase(i,j,2)
+        read(isidep_scbase,*) &
+       dhead_scbasei(i,j), &
+       dhead_scbasej(i,j), &
+       rborn_scbasei(i,j),rborn_scbasej(i,j)
+        read(isidep_scbase,*) &
+       (wdipdip_scbase(k,i,j),k=1,3), &
+       (wqdip_scbase(k,i,j),k=1,2)
+        read(isidep_scbase,*) &
+       alphapol_scbase(i,j), &
+       epsintab_scbase(i,j) 
+       END DO
+      END DO
+      allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
+      allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
+
+      do i=1,ntyp_molec(1)
+       do j=1,ntyp_molec(2)-1 
+          epsij=eps_scbase(i,j)
+          rrij=sigma_scbase(i,j)
+!          r0(i,j)=rrij
+!          r0(j,i)=rrij
+          rrij=rrij**expon
+!          epsij=eps(i,j)
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_scbase(i,j)=epsij*rrij*rrij
+          bb_scbase(i,j)=-sigeps*epsij*rrij
+        enddo
+       enddo
+!-----------------READING PEP BASE POTENTIALS-------------------
+      allocate(eps_pepbase(ntyp_molec(2)))
+      allocate(sigma_pepbase(ntyp_molec(2)))
+      allocate(chi_pepbase(ntyp_molec(2),2))
+      allocate(chipp_pepbase(ntyp_molec(2),2))
+      allocate(alphasur_pepbase(4,ntyp_molec(2)))
+      allocate(sigmap1_pepbase(ntyp_molec(2)))
+      allocate(sigmap2_pepbase(ntyp_molec(2)))
+      allocate(chis_pepbase(ntyp_molec(2),2))
+      allocate(wdipdip_pepbase(3,ntyp_molec(2)))
+
+
+       do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
+        write (*,*) "Im in ", i, " ", j
+        read(isidep_pepbase,*) &
+        eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
+        chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
+         write(*,*) "eps",eps_pepbase(j)
+        read(isidep_pepbase,*) &
+       (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
+       chis_pepbase(j,1),chis_pepbase(j,2)
+        read(isidep_pepbase,*) &
+       (wdipdip_pepbase(k,j),k=1,3)
+       END DO
+      allocate(aa_pepbase(ntyp_molec(2)))
+      allocate(bb_pepbase(ntyp_molec(2)))
+
+       do j=1,ntyp_molec(2)-1
+          epsij=eps_pepbase(j)
+          rrij=sigma_pepbase(j)
+!          r0(i,j)=rrij
+!          r0(j,i)=rrij
+          rrij=rrij**expon
+!          epsij=eps(i,j)
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_pepbase(j)=epsij*rrij*rrij
+          bb_pepbase(j)=-sigeps*epsij*rrij
+        enddo
+!--------------READING SC PHOSPHATE------------------------------------- 
+      allocate(eps_scpho(ntyp_molec(1)))
+      allocate(sigma_scpho(ntyp_molec(1)))
+      allocate(chi_scpho(ntyp_molec(1),2))
+      allocate(chipp_scpho(ntyp_molec(1),2))
+      allocate(alphasur_scpho(4,ntyp_molec(1)))
+      allocate(sigmap1_scpho(ntyp_molec(1)))
+      allocate(sigmap2_scpho(ntyp_molec(1)))
+      allocate(chis_scpho(ntyp_molec(1),2))
+      allocate(wqq_scpho(ntyp_molec(1)))
+      allocate(wqdip_scpho(2,ntyp_molec(1)))
+      allocate(alphapol_scpho(ntyp_molec(1)))
+      allocate(epsintab_scpho(ntyp_molec(1)))
+      allocate(dhead_scphoi(ntyp_molec(1)))
+      allocate(rborn_scphoi(ntyp_molec(1)))
+      allocate(rborn_scphoj(ntyp_molec(1)))
+      allocate(alphi_scpho(ntyp_molec(1)))
+
+
+!      j=1
+       do j=1,ntyp_molec(1) ! without U then we will take T for U
+        write (*,*) "Im in scpho ", i, " ", j
+        read(isidep_scpho,*) &
+        eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
+        chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
+         write(*,*) "eps",eps_scpho(j)
+        read(isidep_scpho,*) &
+       (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
+       chis_scpho(j,1),chis_scpho(j,2)
+        read(isidep_scpho,*) &
+       (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
+        read(isidep_scpho,*) &
+         epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
+         alphi_scpho(j)
+       
+       END DO
+      allocate(aa_scpho(ntyp_molec(1)))
+      allocate(bb_scpho(ntyp_molec(1)))
+
+       do j=1,ntyp_molec(1)
+          epsij=eps_scpho(j)
+          rrij=sigma_scpho(j)
+!          r0(i,j)=rrij
+!          r0(j,i)=rrij
+          rrij=rrij**expon
+!          epsij=eps(i,j)
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_scpho(j)=epsij*rrij*rrij
+          bb_scpho(j)=-sigeps*epsij*rrij
+        enddo
+
+
+        read(isidep_peppho,*) &
+        eps_peppho,sigma_peppho
+        read(isidep_peppho,*) &
+       (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
+        read(isidep_peppho,*) &
+       (wqdip_peppho(k),k=1,2)
+
+          epsij=eps_peppho
+          rrij=sigma_peppho
+!          r0(i,j)=rrij
+!          r0(j,i)=rrij
+          rrij=rrij**expon
+!          epsij=eps(i,j)
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_peppho=epsij*rrij*rrij
+          bb_peppho=-sigeps*epsij*rrij
+

 
       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
       bad(:,:)=0.0D0
 
       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
       bad(:,:)=0.0D0

@@ -2216,6 +2721,20 @@
       endif
 !      lprint=.false.
 #endif
       endif
 !      lprint=.false.
 #endif

+      allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
+
+      do i=1,ntyp_molec(2)
+        read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
+      enddo
+      do i=1,ntyp_molec(2)
+        aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
+        bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
+      enddo
+      r0pp=1.12246204830937298142*5.16158
+      epspp=4.95713/4
+      AEES=108.661
+      BEES=0.433246
+

 !
 ! Define the constants of the disulfide bridge
 !
 !
 ! Define the constants of the disulfide bridge
 !

@@ -2252,6 +2771,31 @@
       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
         ' v3ss:',v3ss
       endif
       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
         ' v3ss:',v3ss
       endif

+      if (shield_mode.gt.0) then
+      pi=4.0D0*datan(1.0D0)
+!C VSolvSphere the volume of solving sphere
+      print *,pi,"pi"
+!C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
+!C there will be no distinction between proline peptide group and normal peptide
+!C group in case of shielding parameters
+      VSolvSphere=4.0/3.0*pi*(4.50d0)**3
+      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
+      write (iout,*) VSolvSphere,VSolvSphere_div
+!C long axis of side chain 
+      do i=1,ntyp
+      long_r_sidechain(i)=vbldsc0(1,i)
+!      if (scelemode.eq.0) then
+      short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
+      if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
+!      else
+!      short_r_sidechain(i)=sigma(i,i)
+!      endif
+      write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
+         sigma0(i) 
+      enddo
+      buff_shield=1.0d0
+      endif
+

       return
   111 write (iout,*) "Error reading bending energy parameters."
       goto 999
       return
   111 write (iout,*) "Error reading bending energy parameters."
       goto 999

@@ -2310,11 +2854,11 @@
 ! Read the PDB file and convert the peptide geometry into virtual-chain 
 ! geometry.
       use geometry_data
 ! Read the PDB file and convert the peptide geometry into virtual-chain 
 ! geometry.
       use geometry_data

-      use energy_data, only: itype
+      use energy_data, only: itype,istype

       use control_data
       use compare_data
       use MPI_data
       use control_data
       use compare_data
       use MPI_data

-      use control, only: rescode
+!      use control, only: rescode,sugarcode

 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.LOCAL'
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.LOCAL'

@@ -2327,7 +2871,7 @@
 !      include 'COMMON.CONTROL'
 !      include 'COMMON.DISTFIT'
 !      include 'COMMON.SETUP'
 !      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,k!,ity!,&

 !        ishift_pdb
       logical :: lprn=.true.,fail
       real(kind=8),dimension(3) :: e1,e2,e3
 !        ishift_pdb
       logical :: lprn=.true.,fail
       real(kind=8),dimension(3) :: e1,e2,e3

@@ -2336,32 +2880,51 @@
       character(len=5) :: atom
       character(len=80) :: card
       real(kind=8),dimension(3,20) :: sccor
       character(len=5) :: atom
       character(len=80) :: card
       real(kind=8),dimension(3,20) :: sccor

-      integer :: kkk,lll,icha,kupa     !rescode,
+      integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin      !rescode,
+      integer :: isugar,molecprev,firstion
+      character*1 :: sugar

       real(kind=8) :: cou
       real(kind=8) :: cou

+      real(kind=8),dimension(3) :: ccc

 !el local varables
       integer,dimension(2,maxres/3) :: hfrag_alloc
       integer,dimension(4,maxres/3) :: bfrag_alloc
       real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
 !el local varables
       integer,dimension(2,maxres/3) :: hfrag_alloc
       integer,dimension(4,maxres/3) :: bfrag_alloc
       real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)

-
+      real(kind=8),dimension(:,:), allocatable  :: c_temporary
+      integer,dimension(:,:) , allocatable  :: itype_temporary
+      integer,dimension(:),allocatable :: istype_temp

       efree_temp=0.0d0
       ibeg=1
       ishift1=0
       ishift=0
       efree_temp=0.0d0
       ibeg=1
       ishift1=0
       ishift=0

+      molecule=1
+      counter=0

 !      write (2,*) "UNRES_PDB",unres_pdb
       ires=0
       ires_old=0
 !      write (2,*) "UNRES_PDB",unres_pdb
       ires=0
       ires_old=0

+#ifdef WHAM_RUN
+      do i=1,nres
+       do j=1,5
+        itype(i,j)=0
+       enddo
+      enddo
+#endif

       nres=0
       iii=0
       lsecondary=.false.
       nhfrag=0
       nbfrag=0
       nres=0
       iii=0
       lsecondary=.false.
       nhfrag=0
       nbfrag=0

+      do j=1,5
+       nres_molec(j)=0
+      enddo
+      
+       

 !-----------------------------
       allocate(hfrag(2,maxres/3)) !(2,maxres/3)
       allocate(bfrag(4,maxres/3)) !(4,maxres/3)
 !-----------------------------
       allocate(hfrag(2,maxres/3)) !(2,maxres/3)
       allocate(bfrag(4,maxres/3)) !(4,maxres/3)

-
+      if(.not. allocated(istype)) allocate(istype(maxres))

       do i=1,100000
         read (ipdbin,'(a80)',end=10) card
       do i=1,100000
         read (ipdbin,'(a80)',end=10) card

-!       write (iout,'(a)') card
+       write (iout,'(a)') card

         if (card(:5).eq.'HELIX') then
           nhfrag=nhfrag+1
           lsecondary=.true.
         if (card(:5).eq.'HELIX') then
           nhfrag=nhfrag+1
           lsecondary=.true.

@@ -2384,9 +2947,11 @@
         else if (card(:3).eq.'TER') then
 ! End current chain
           ires_old=ires+2
         else if (card(:3).eq.'TER') then
 ! End current chain
           ires_old=ires+2

+          ishift=ishift+1

           ishift1=ishift1+1
           ishift1=ishift1+1

-          itype(ires_old)=ntyp1
-          itype(ires_old-1)=ntyp1
+          itype(ires_old,molecule)=ntyp1_molec(molecule)
+          itype(ires_old-1,molecule)=ntyp1_molec(molecule)
+          nres_molec(molecule)=nres_molec(molecule)+2

           ibeg=2
 !          write (iout,*) "Chain ended",ires,ishift,ires_old
           if (unres_pdb) then
           ibeg=2
 !          write (iout,*) "Chain ended",ires,ishift,ires_old
           if (unres_pdb) then

@@ -2402,6 +2967,7 @@
 ! Read free energy
         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
 ! Fish out the ATOM cards.
 ! Read free energy
         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
 ! Fish out the ATOM cards.

+!        write(iout,*) 'card',card(1:20)

         if (index(card(1:4),'ATOM').gt.0) then  
           read (card(12:16),*) atom
 !          write (iout,*) "! ",atom," !",ires
         if (index(card(1:4),'ATOM').gt.0) then  
           read (card(12:16),*) atom
 !          write (iout,*) "! ",atom," !",ires

@@ -2435,7 +3001,8 @@
               ishift=ires-1
               if (res.ne.'GLY' .and. res.ne. 'ACE') then
                 ishift=ishift-1
               ishift=ires-1
               if (res.ne.'GLY' .and. res.ne. 'ACE') then
                 ishift=ishift-1

-                itype(1)=ntyp1
+                itype(1,1)=ntyp1
+                nres_molec(molecule)=nres_molec(molecule)+1

               endif
               ires=ires-ishift+ishift1
               ires_old=ires
               endif
               ires=ires-ishift+ishift1
               ires_old=ires

@@ -2448,17 +3015,37 @@
               ishift1=ishift1-1    !!!!!
 !              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
               ires=ires-ishift+ishift1
               ishift1=ishift1-1    !!!!!
 !              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
               ires=ires-ishift+ishift1

+!              print *,ires,ishift,ishift1

               ires_old=ires
               ibeg=0
             else
               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
               ires=ires-ishift+ishift1
               ires_old=ires
               ires_old=ires
               ibeg=0
             else
               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
               ires=ires-ishift+ishift1
               ires_old=ires

-            endif
+            endif 
+!            print *,'atom',ires,atom

             if (res.eq.'ACE' .or. res.eq.'NHE') then
             if (res.eq.'ACE' .or. res.eq.'NHE') then

-              itype(ires)=10
+              itype(ires,1)=10
+            else
+             if (atom.eq.'CA  '.or.atom.eq.'N   ') then
+             molecule=1
+              itype(ires,molecule)=rescode(ires,res,0,molecule)
+              firstion=0
+!              nres_molec(molecule)=nres_molec(molecule)+1

             else
             else

-              itype(ires)=rescode(ires,res,0)
+             molecule=2
+              itype(ires,molecule)=rescode(ires,res(2:3),0,molecule)
+!              nres_molec(molecule)=nres_molec(molecule)+1
+             read (card(19:19),'(a1)') sugar
+             isugar=sugarcode(sugar,ires)
+!            if (ibeg.eq.1) then
+!              istype(1)=isugar
+!            else
+              istype(ires)=isugar
+!              print *,"ires=",ires,istype(ires)
+!            endif
+
+            endif

             endif
           else
             ires=ires-ishift+ishift1
             endif
           else
             ires=ires-ishift+ishift1

@@ -2471,91 +3058,197 @@
           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)
           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)

+!              print *,ires,ishift,ishift1

 !            write (iout,*) "backbone ",atom
 #ifdef DEBUG
             write (iout,'(2i3,2x,a,3f8.3)') &
 !            write (iout,*) "backbone ",atom
 #ifdef DEBUG
             write (iout,'(2i3,2x,a,3f8.3)') &

-            ires,itype(ires),res,(c(j,ires),j=1,3)
+            ires,itype(ires,1),res,(c(j,ires),j=1,3)

 #endif
             iii=iii+1
 #endif
             iii=iii+1

+              nres_molec(molecule)=nres_molec(molecule)+1

             do j=1,3
               sccor(j,iii)=c(j,ires)
             enddo
             do j=1,3
               sccor(j,iii)=c(j,ires)
             enddo

-!            write (*,*) card(23:27),ires,itype(ires)
+          else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
+               atom.eq."C2'" .or. atom.eq."C3'" &
+               .or. atom.eq."C4'" .or. atom.eq."O4'")) then
+            read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
+!c            write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
+!              print *,ires,ishift,ishift1
+            counter=counter+1
+!            iii=iii+1
+!            do j=1,3
+!              sccor(j,iii)=c(j,ires)
+!            enddo
+            do j=1,3
+              c(j,ires)=c(j,ires)+ccc(j)/5.0
+            enddo
+             print *,counter,molecule
+             if (counter.eq.5) then
+!            iii=iii+1
+              nres_molec(molecule)=nres_molec(molecule)+1
+              firstion=0
+!            do j=1,3
+!              sccor(j,iii)=c(j,ires)
+!            enddo
+             counter=0
+           endif
+!            print *, "ATOM",atom(1:3)
+          else if (atom.eq."C5'") then
+             read (card(19:19),'(a1)') sugar
+             isugar=sugarcode(sugar,ires)
+            if (ibeg.eq.1) then
+              istype(1)=isugar
+            else
+              istype(ires)=isugar
+!              print *,ires,istype(ires)
+            endif
+            if (unres_pdb) then
+              molecule=2
+!              print *,"nres_molec(molecule)",nres_molec(molecule),ires
+              read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+              nres_molec(molecule)=nres_molec(molecule)+1
+              print *,"nres_molec(molecule)",nres_molec(molecule),ires
+
+            else
+              iii=iii+1
+              read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+            endif
+          else if ((atom.eq."C1'").and.unres_pdb) then
+              iii=iii+1
+              read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+!            write (*,*) card(23:27),ires,itype(ires,1)

           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
                    atom.ne.'N' .and. atom.ne.'C' .and. &
                    atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
                    atom.ne.'N' .and. atom.ne.'C' .and. &
                    atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &

-                   atom.ne.'OXT' .and. atom(:2).ne.'3H') then
+                   atom.ne.'OXT' .and. atom(:2).ne.'3H' &
+                   .and. atom.ne.'P  '.and. &
+                  atom(1:1).ne.'H' .and. &
+                  atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
+                  ) then

 !            write (iout,*) "sidechain ",atom
 !            write (iout,*) "sidechain ",atom

+!            write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
+                 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
+!                        write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)

             iii=iii+1
             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
             iii=iii+1
             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)

+              endif

           endif
           endif

-        endif
+        else if ((ions).and.(card(1:6).eq.'HETATM')) then
+       if (firstion.eq.0) then 
+       firstion=1
+       if (unres_pdb) then
+         do j=1,3
+           dc(j,ires)=sccor(j,iii)
+         enddo
+       else
+          call sccenter(ires,iii,sccor)
+       endif
+       endif
+          read (card(12:16),*) atom
+          print *,"HETATOM", atom
+          read (card(18:20),'(a3)') res
+          if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
+          (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG')           &
+          .or.(atom(1:2).eq.'K ')) &
+          then
+           ires=ires+1
+           if (molecule.ne.5) molecprev=molecule
+           molecule=5
+           nres_molec(molecule)=nres_molec(molecule)+1
+           print *,"HERE",nres_molec(molecule)
+           res=res(2:3)//' '
+           itype(ires,molecule)=rescode(ires,res,0,molecule)
+           read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+          endif
+        endif !atom

       enddo
    10 write (iout,'(a,i5)') ' Number of residues found: ',ires
       if (ires.eq.0) return
 ! Calculate dummy residue coordinates inside the "chain" of a multichain
 ! system
       nres=ires
       enddo
    10 write (iout,'(a,i5)') ' Number of residues found: ',ires
       if (ires.eq.0) return
 ! Calculate dummy residue coordinates inside the "chain" of a multichain
 ! system
       nres=ires

+      if (((ires_old.ne.ires).and.(molecule.ne.5)) &
+        ) &
+         nres_molec(molecule)=nres_molec(molecule)-2
+      print *,'I have',nres, nres_molec(:)
+      
+      do k=1,4 ! ions are without dummy 
+       if (nres_molec(k).eq.0) cycle

       do i=2,nres-1
       do i=2,nres-1

-!        write (iout,*) i,itype(i)
-!        if (itype(i).eq.ntyp1) then
-!          write (iout,*) "dummy",i,itype(i)
+!        write (iout,*) i,itype(i,1)
+!        if (itype(i,1).eq.ntyp1) then
+!          write (iout,*) "dummy",i,itype(i,1)

 !          do j=1,3
 !            c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
 !            c(j,i)=(c(j,i-1)+c(j,i+1))/2
 !            dc(j,i)=c(j,i)
 !          enddo
 !        endif
 !          do j=1,3
 !            c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
 !            c(j,i)=(c(j,i-1)+c(j,i+1))/2
 !            dc(j,i)=c(j,i)
 !          enddo
 !        endif

-        if (itype(i).eq.ntyp1) then
-         if (itype(i+1).eq.ntyp1) then
+        if (itype(i,k).eq.ntyp1_molec(k)) then
+         if (itype(i+1,k).eq.ntyp1_molec(k)) then
+          if (itype(i+2,k).eq.0) then 
+!           print *,"masz sieczke"
+           do j=1,5
+            if (itype(i+2,j).ne.0) then
+            itype(i+1,k)=0
+            itype(i+1,j)=ntyp1_molec(j)
+            nres_molec(k)=nres_molec(k)-1
+            nres_molec(j)=nres_molec(j)+1
+            go to 3331
+            endif
+           enddo
+ 3331      continue
+          endif

 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain

-! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
-! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
-           if (unres_pdb) then
+! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
+! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
+!           if (unres_pdb) then

 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
 !            print *,i,'tu dochodze'
 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
 !            print *,i,'tu dochodze'

-            call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
-            if (fail) then
-              e2(1)=0.0d0
-              e2(2)=1.0d0
-              e2(3)=0.0d0
-            endif !fail
-            print *,i,'a tu?'
-            do j=1,3
-             c(j,i)=c(j,i-1)-1.9d0*e2(j)
-            enddo
-           else   !unres_pdb
+!            call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
+!            if (fail) then
+!              e2(1)=0.0d0
+!              e2(2)=1.0d0
+!              e2(3)=0.0d0
+!            endif !fail
+!            print *,i,'a tu?'
+!            do j=1,3
+!             c(j,i)=c(j,i-1)-1.9d0*e2(j)
+!            enddo
+!           else   !unres_pdb

            do j=1,3
              dcj=(c(j,i-2)-c(j,i-3))/2.0
             if (dcj.eq.0) dcj=1.23591524223
              c(j,i)=c(j,i-1)+dcj
              c(j,nres+i)=c(j,i)
            enddo
            do j=1,3
              dcj=(c(j,i-2)-c(j,i-3))/2.0
             if (dcj.eq.0) dcj=1.23591524223
              c(j,i)=c(j,i-1)+dcj
              c(j,nres+i)=c(j,i)
            enddo

-          endif   !unres_pdb
-         else     !itype(i+1).eq.ntyp1
-          if (unres_pdb) then
+!          endif   !unres_pdb
+         else     !itype(i+1,1).eq.ntyp1
+!          if (unres_pdb) then

 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue

-            call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
-            if (fail) then
-              e2(1)=0.0d0
-              e2(2)=1.0d0
-              e2(3)=0.0d0
-            endif
+!            call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
+!            if (fail) then
+!              e2(1)=0.0d0
+!              e2(2)=1.0d0
+!              e2(3)=0.0d0
+!            endif

             do j=1,3
               c(j,i)=c(j,i+1)-1.9d0*e2(j)
             enddo
             do j=1,3
               c(j,i)=c(j,i+1)-1.9d0*e2(j)
             enddo

-          else !unres_pdb
+!          else !unres_pdb

            do j=1,3
             dcj=(c(j,i+3)-c(j,i+2))/2.0
             if (dcj.eq.0) dcj=1.23591524223
             c(j,i)=c(j,i+1)-dcj
             c(j,nres+i)=c(j,i)
            enddo
            do j=1,3
             dcj=(c(j,i+3)-c(j,i+2))/2.0
             if (dcj.eq.0) dcj=1.23591524223
             c(j,i)=c(j,i+1)-dcj
             c(j,nres+i)=c(j,i)
            enddo

-          endif !unres_pdb
-         endif !itype(i+1).eq.ntyp1
+!          endif !unres_pdb
+         endif !itype(i+1,1).eq.ntyp1

         endif  !itype.eq.ntyp1
 
       enddo
         endif  !itype.eq.ntyp1
 
       enddo

+      enddo

 ! Calculate the CM of the last side chain.
       if (iii.gt.0)  then
       if (unres_pdb) then
 ! Calculate the CM of the last side chain.
       if (iii.gt.0)  then
       if (unres_pdb) then

@@ -2569,28 +3262,33 @@
 !      nres=ires
       nsup=nres
       nstart_sup=1
 !      nres=ires
       nsup=nres
       nstart_sup=1

-      if (itype(nres).ne.10) then
+!      print *,"molecule",molecule
+      if ((itype(nres,1).ne.10)) then

         nres=nres+1
         nres=nres+1

-        itype(nres)=ntyp1
-        if (unres_pdb) then
+          if (molecule.eq.5) molecule=molecprev
+        itype(nres,molecule)=ntyp1_molec(molecule)
+        nres_molec(molecule)=nres_molec(molecule)+1
+!        if (unres_pdb) then

 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue

-          call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
-          if (fail) then
-            e2(1)=0.0d0
-            e2(2)=1.0d0
-            e2(3)=0.0d0
-          endif
-          do j=1,3
-            c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
-          enddo
-        else
+!          call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
+!          if (fail) then
+!            e2(1)=0.0d0
+!            e2(2)=1.0d0
+!            e2(3)=0.0d0
+!          endif
+!          do j=1,3
+!            c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
+!          enddo
+!        else

         do j=1,3
           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
           c(j,nres)=c(j,nres-1)+dcj
           c(j,2*nres)=c(j,nres)
         enddo
         do j=1,3
           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
           c(j,nres)=c(j,nres-1)+dcj
           c(j,2*nres)=c(j,nres)
         enddo

-        endif
+!        endif

       endif
       endif

+!     print *,'I have',nres, nres_molec(:)
+

 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
 #ifdef WHAM_RUN
       if (nres.ne.nres0) then
 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
 #ifdef WHAM_RUN
       if (nres.ne.nres0) then

@@ -2638,7 +3336,8 @@
         c(j,nres+1)=c(j,1)
         c(j,2*nres)=c(j,nres)
       enddo
         c(j,nres+1)=c(j,1)
         c(j,2*nres)=c(j,nres)
       enddo

-      if (itype(1).eq.ntyp1) then
+      
+      if (itype(1,1).eq.ntyp1) then

         nsup=nsup-1
         nstart_sup=2
         if (unres_pdb) then
         nsup=nsup-1
         nstart_sup=2
         if (unres_pdb) then

@@ -2660,6 +3359,24 @@
         enddo
         endif
       endif
         enddo
         endif
       endif

+! First lets assign correct dummy to correct type of chain
+! 1) First residue
+      if (itype(1,1).eq.ntyp1) then
+        if (itype(2,1).eq.0) then
+         do j=2,5
+           if (itype(2,j).ne.0) then
+           itype(1,1)=0
+           itype(1,j)=ntyp1_molec(j)
+           nres_molec(1)=nres_molec(1)-1
+           nres_molec(j)=nres_molec(j)+1
+           go to 3231
+           endif
+         enddo
+3231    continue
+        endif
+       endif
+       print *,'I have',nres, nres_molec(:)
+

 ! Copy the coordinates to reference coordinates
 !      do i=1,2*nres
 !        do j=1,3
 ! Copy the coordinates to reference coordinates
 !      do i=1,2*nres
 !        do j=1,3

@@ -2670,11 +3387,11 @@
       if (lprn) then
       write (iout,'(/a)') &
         "Cartesian coordinates of the reference structure"
       if (lprn) then
       write (iout,'(/a)') &
         "Cartesian coordinates of the reference structure"

-      write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
+      write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &

        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
       do ires=1,nres
        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
       do ires=1,nres

-        write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
-          restyp(itype(ires)),ires,(c(j,ires),j=1,3),&
+        write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
+          (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&

           (c(j,ires+nres),j=1,3)
       enddo
       endif
           (c(j,ires+nres),j=1,3)
       enddo
       endif

@@ -2684,12 +3401,83 @@
        write (iout,'(a)') &
          "Backbone and SC coordinates as read from the PDB"
        do ires=1,nres
        write (iout,'(a)') &
          "Backbone and SC coordinates as read from the PDB"
        do ires=1,nres

-        write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
-          ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),&
+        write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
+          ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&

           (c(j,nres+ires),j=1,3)
        enddo
       endif
           (c(j,nres+ires),j=1,3)
        enddo
       endif

+! NOW LETS ROCK! SORTING
+      allocate(c_temporary(3,2*nres))
+      allocate(itype_temporary(nres,5))
+      if (.not.allocated(molnum)) allocate(molnum(nres+1))
+      if (.not.allocated(istype)) write(iout,*) &
+          "SOMETHING WRONG WITH ISTYTPE"
+      allocate(istype_temp(nres))
+       itype_temporary(:,:)=0
+      seqalingbegin=1
+      do k=1,5
+        do i=1,nres
+         if (itype(i,k).ne.0) then
+          do j=1,3
+          c_temporary(j,seqalingbegin)=c(j,i)
+          c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
+
+          enddo
+          itype_temporary(seqalingbegin,k)=itype(i,k)
+          print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
+          istype_temp(seqalingbegin)=istype(i)
+          molnum(seqalingbegin)=k
+          seqalingbegin=seqalingbegin+1
+         endif
+        enddo
+       enddo
+       do i=1,2*nres
+        do j=1,3
+        c(j,i)=c_temporary(j,i)
+        enddo
+       enddo
+       do k=1,5
+        do i=1,nres
+         itype(i,k)=itype_temporary(i,k)
+         istype(i)=istype_temp(i)
+        enddo
+       enddo
+!      if (itype(1,1).eq.ntyp1) then
+!        nsup=nsup-1
+!        nstart_sup=2
+!        if (unres_pdb) then
+! 2/15/2013 by Adam: corrected insertion of the first dummy residue
+!          call refsys(2,3,4,e1,e2,e3,fail)
+!          if (fail) then
+!            e2(1)=0.0d0
+!            e2(2)=1.0d0
+!            e2(3)=0.0d0
+!          endif
+!          do j=1,3
+!            c(j,1)=c(j,2)-1.9d0*e2(j)
+!          enddo
+!        else
+!        do j=1,3
+!          dcj=(c(j,4)-c(j,3))/2.0
+!          c(j,1)=c(j,2)-dcj
+!          c(j,nres+1)=c(j,1)
+!        enddo
+!        endif
+!      endif
+
+      if (lprn) then
+      write (iout,'(/a)') &
+        "Cartesian coordinates of the reference structure after sorting"
+      write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
+       "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+      do ires=1,nres
+        write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
+          (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
+          (c(j,ires+nres),j=1,3)
+      enddo
+      endif

 
 

+!       print *,seqalingbegin,nres

       if(.not.allocated(vbld)) then
        allocate(vbld(2*nres))
        do i=1,2*nres
       if(.not.allocated(vbld)) then
        allocate(vbld(2*nres))
        do i=1,2*nres

@@ -2768,11 +3556,12 @@
       kkk=1
       lll=0
       cou=1
       kkk=1
       lll=0
       cou=1

+        write (iout,*) "symetr", symetr

       do i=1,nres
       lll=lll+1
       do i=1,nres
       lll=lll+1

-!c      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
+!      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)

       if (i.gt.1) then
       if (i.gt.1) then

-      if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
+      if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then

       chain_length=lll-1
       kkk=kkk+1
 !       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
       chain_length=lll-1
       kkk=kkk+1
 !       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)

@@ -2799,7 +3588,7 @@
 !       write (iout,*) "spraw lancuchy",chain_length,symetr
 !       do i=1,4
 !         do kkk=1,chain_length
 !       write (iout,*) "spraw lancuchy",chain_length,symetr
 !       do i=1,4
 !         do kkk=1,chain_length

-!           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
+!           write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)

 !         enddo
 !        enddo
 ! enddiagnostic
 !         enddo
 !        enddo
 ! enddiagnostic

@@ -2822,7 +3611,7 @@
         cou=0
         do kkk=1,symetr
          icha=tabperm(i,kkk)
         cou=0
         do kkk=1,symetr
          icha=tabperm(i,kkk)

-!         write (iout,*) i,icha
+         write (iout,*) i,icha

          do lll=1,chain_length
           cou=cou+1
            if (cou.le.nres) then
          do lll=1,chain_length
           cou=cou+1
            if (cou.le.nres) then

@@ -2830,7 +3619,7 @@
             kupa=mod(lll,chain_length)
             iprzes=(kkk-1)*chain_length+lll
             if (kupa.eq.0) kupa=chain_length
             kupa=mod(lll,chain_length)
             iprzes=(kkk-1)*chain_length+lll
             if (kupa.eq.0) kupa=chain_length

-!            write (iout,*) "kupa", kupa
+            write (iout,*) "kupa", kupa

             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
           enddo
             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
           enddo

@@ -2844,16 +3633,16 @@
       do kkk=1,nperm
       write (iout,*) "nowa struktura", nperm
       do i=1,nres
       do kkk=1,nperm
       write (iout,*) "nowa struktura", nperm
       do i=1,nres

-        write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),&
+        write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&

       cref(2,i,kkk),&
       cref(3,i,kkk),cref(1,nres+i,kkk),&
       cref(2,nres+i,kkk),cref(3,nres+i,kkk)
       enddo
       cref(2,i,kkk),&
       cref(3,i,kkk),cref(1,nres+i,kkk),&
       cref(2,nres+i,kkk),cref(3,nres+i,kkk)
       enddo

-  100 format (//'              alpha-carbon coordinates       ',&
+  100 format (//'                alpha-carbon coordinates       ',&

                 '     centroid coordinates'/ &
                 '       ', 6X,'X',11X,'Y',11X,'Z', &
                                 10X,'X',11X,'Y',11X,'Z')
                 '     centroid coordinates'/ &
                 '       ', 6X,'X',11X,'Y',11X,'Z', &
                                 10X,'X',11X,'Y',11X,'Z')

-  110 format (a,'(',i3,')',6f12.5)
+  110 format (a,'(',i5,')',6f12.5)

      
       enddo
 !c enddiag
      
       enddo
 !c enddiag

@@ -2970,10 +3759,26 @@
 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
         write(iout,*) "shield_mode",shield_mode
 !C  Varibles set size of box
 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
         write(iout,*) "shield_mode",shield_mode
 !C  Varibles set size of box

+      with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
+      protein=index(controlcard,"PROTEIN").gt.0
+      ions=index(controlcard,"IONS").gt.0
+      nucleic=index(controlcard,"NUCLEIC").gt.0
+      write (iout,*) "with_theta_constr ",with_theta_constr
+      AFMlog=(index(controlcard,'AFM'))
+      selfguide=(index(controlcard,'SELFGUIDE'))
+      print *,'AFMlog',AFMlog,selfguide,"KUPA"
+      call readi(controlcard,'GENCONSTR',genconstr,0)

       call reada(controlcard,'BOXX',boxxsize,100.0d0)
       call reada(controlcard,'BOXY',boxysize,100.0d0)
       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
       call readi(controlcard,'TUBEMOD',tubemode,0)
       call reada(controlcard,'BOXX',boxxsize,100.0d0)
       call reada(controlcard,'BOXY',boxysize,100.0d0)
       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
       call readi(controlcard,'TUBEMOD',tubemode,0)

+      print *,"SCELE",scelemode
+      call readi(controlcard,"SCELEMODE",scelemode,0)
+      print *,"SCELE",scelemode
+
+! elemode = 0 is orignal UNRES electrostatics
+! elemode = 1 is "Momo" potentials in progress
+! elemode = 2 is in development EVALD

       write (iout,*) TUBEmode,"TUBEMODE"
       if (TUBEmode.gt.0) then
        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
       write (iout,*) TUBEmode,"TUBEMODE"
       if (TUBEmode.gt.0) then
        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)

@@ -3111,22 +3916,22 @@
        write (iout,'(2a)') diagmeth(kdiag),&
         ' routine used to diagonalize matrices.'
       if (shield_mode.gt.0) then
        write (iout,'(2a)') diagmeth(kdiag),&
         ' routine used to diagonalize matrices.'
       if (shield_mode.gt.0) then

-      pi=3.141592d0
+       pi=4.0D0*datan(1.0D0)

 !C VSolvSphere the volume of solving sphere
 !C VSolvSphere the volume of solving sphere

-!C      print *,pi,"pi"
+      print *,pi,"pi"

 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
 !C there will be no distinction between proline peptide group and normal peptide
 !C group in case of shielding parameters
 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
 !C there will be no distinction between proline peptide group and normal peptide
 !C group in case of shielding parameters

-      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
-      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
+      VSolvSphere=4.0/3.0*pi*(4.50d0)**3
+      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3

       write (iout,*) VSolvSphere,VSolvSphere_div
 !C long axis of side chain 
       write (iout,*) VSolvSphere,VSolvSphere_div
 !C long axis of side chain 

-      do i=1,ntyp
-      long_r_sidechain(i)=vbldsc0(1,i)
-      short_r_sidechain(i)=sigma0(i)
-      write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
-         sigma0(i) 
-      enddo
+!      do i=1,ntyp
+!      long_r_sidechain(i)=vbldsc0(1,i)
+!      short_r_sidechain(i)=sigma0(i)
+!      write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
+!         sigma0(i) 
+!      enddo

       buff_shield=1.0d0
       endif
       return
       buff_shield=1.0d0
       endif
       return

@@ -3251,7 +4056,7 @@
 !      character(len=80) :: ucase
       character(len=320) :: controlcard
 !el local variables
 !      character(len=80) :: ucase
       character(len=320) :: controlcard
 !el local variables

-      integer :: i
+      integer :: i,j

       real(kind=8) :: eta
 
       call card_concat(controlcard,.true.)
       real(kind=8) :: eta
 
       call card_concat(controlcard,.true.)

@@ -3286,6 +4091,10 @@
       large = index(controlcard,"LARGE").gt.0
       print_compon = index(controlcard,"PRINT_COMPON").gt.0
       rattle = index(controlcard,"RATTLE").gt.0
       large = index(controlcard,"LARGE").gt.0
       print_compon = index(controlcard,"PRINT_COMPON").gt.0
       rattle = index(controlcard,"RATTLE").gt.0

+      preminim=(index(controlcard,'PREMINIM').gt.0)
+      write (iout,*) "PREMINIM ",preminim
+      dccart=(index(controlcard,'CART').gt.0)
+      if (preminim) call read_minim

 !  if performing umbrella sampling, fragments constrained are read from the fragment file 
       nset=0
       if(usampl) then
 !  if performing umbrella sampling, fragments constrained are read from the fragment file 
       nset=0
       if(usampl) then

@@ -3364,21 +4173,27 @@
         if(me.eq.king.or..not.out1file) &
          write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
           eta
         if(me.eq.king.or..not.out1file) &
          write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
           eta

-        gamp=scal_fric*(pstok+rwat)*eta
-        stdfp=dsqrt(2*Rb*t_bath/d_time)
-        allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
+!        allocate(gamp
+        do j=1,5 !types of molecules
+        gamp(j)=scal_fric*(pstok(j)+rwat)*eta
+        stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
+        enddo
+        allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
+        do j=1,5 !types of molecules

         do i=1,ntyp
         do i=1,ntyp

-          gamsc(i)=scal_fric*(restok(i)+rwat)*eta  
-          stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
+          gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta  
+          stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)

         enddo 
         enddo 

+        enddo
+

         if(me.eq.king.or..not.out1file)then
          write (iout,'(/2a/)') &
          "Radii of site types and friction coefficients and std's of",&
          " stochastic forces of fully exposed sites"
         if(me.eq.king.or..not.out1file)then
          write (iout,'(/2a/)') &
          "Radii of site types and friction coefficients and std's of",&
          " stochastic forces of fully exposed sites"

-         write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
+         write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))

          do i=1,ntyp
          do i=1,ntyp

-          write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),&
-           gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
+          write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
+           gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))

          enddo
         endif
       else if (tbf) then
          enddo
         endif
       else if (tbf) then

@@ -3739,6 +4554,8 @@
 ! Get parameter filenames and open the parameter files.
       call getenv_loc('BONDPAR',bondname)
       open (ibond,file=bondname,status='old',readonly,shared)
 ! Get parameter filenames and open the parameter files.
       call getenv_loc('BONDPAR',bondname)
       open (ibond,file=bondname,status='old',readonly,shared)

+      call getenv_loc('BONDPAR_NUCL',bondname_nucl)
+      open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)

       call getenv_loc('THETPAR',thetname)
       open (ithep,file=thetname,status='old',readonly,shared)
       call getenv_loc('ROTPAR',rotname)
       call getenv_loc('THETPAR',thetname)
       open (ithep,file=thetname,status='old',readonly,shared)
       call getenv_loc('ROTPAR',rotname)

@@ -3753,6 +4570,19 @@
       open (ielep,file=elename,status='old',readonly,shared)
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old',readonly,shared)
       open (ielep,file=elename,status='old',readonly,shared)
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old',readonly,shared)

+
+      call getenv_loc('THETPAR_NUCL',thetname_nucl)
+      open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
+      call getenv_loc('ROTPAR_NUCL',rotname_nucl)
+      open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
+      call getenv_loc('TORPAR_NUCL',torname_nucl)
+      open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
+      call getenv_loc('TORDPAR_NUCL',tordname_nucl)
+      open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
+      call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
+      open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
+
+

 #elif (defined CRAY) || (defined AIX)
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
         action='read')
 #elif (defined CRAY) || (defined AIX)
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
         action='read')

@@ -3763,6 +4593,9 @@
 ! Get parameter filenames and open the parameter files.
       call getenv_loc('BONDPAR',bondname)
       open (ibond,file=bondname,status='old',action='read')
 ! Get parameter filenames and open the parameter files.
       call getenv_loc('BONDPAR',bondname)
       open (ibond,file=bondname,status='old',action='read')

+      call getenv_loc('BONDPAR_NUCL',bondname_nucl)
+      open (ibond_nucl,file=bondname_nucl,status='old',action='read')
+

 !      print *,"Processor",myrank," opened file IBOND" 
       call getenv_loc('THETPAR',thetname)
       open (ithep,file=thetname,status='old',action='read')
 !      print *,"Processor",myrank," opened file IBOND" 
       call getenv_loc('THETPAR',thetname)
       open (ithep,file=thetname,status='old',action='read')

@@ -3787,10 +4620,24 @@
 !      print *,"Processor",myrank," opened file IELEP" 
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old',action='read')
 !      print *,"Processor",myrank," opened file IELEP" 
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old',action='read')

+
+      call getenv_loc('THETPAR_NUCL',thetname_nucl)
+      open (ithep_nucl,file=thetname_nucl,status='old',action='read')
+      call getenv_loc('ROTPAR_NUCL',rotname_nucl)
+      open (irotam_nucl,file=rotname_nucl,status='old',action='read')
+      call getenv_loc('TORPAR_NUCL',torname_nucl)
+      open (itorp_nucl,file=torname_nucl,status='old',action='read')
+      call getenv_loc('TORDPAR_NUCL',tordname_nucl)
+      open (itordp_nucl,file=tordname_nucl,status='old',action='read')
+      call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
+      open (isidep_nucl,file=sidename_nucl,status='old',action='read')
+

       call getenv_loc('LIPTRANPAR',liptranname)
       open (iliptranpar,file=liptranname,status='old',action='read')
       call getenv_loc('TUBEPAR',tubename)
       open (itube,file=tubename,status='old',action='read')
       call getenv_loc('LIPTRANPAR',liptranname)
       open (iliptranpar,file=liptranname,status='old',action='read')
       call getenv_loc('TUBEPAR',tubename)
       open (itube,file=tubename,status='old',action='read')

+      call getenv_loc('IONPAR',ionname)
+      open (iion,file=ionname,status='old',action='read')

 
 !      print *,"Processor",myrank," opened file ISIDEP" 
 !      print *,"Processor",myrank," opened parameter files" 
 
 !      print *,"Processor",myrank," opened file ISIDEP" 
 !      print *,"Processor",myrank," opened parameter files" 

@@ -3801,6 +4648,9 @@
 ! Get parameter filenames and open the parameter files.
       call getenv_loc('BONDPAR',bondname)
       open (ibond,file=bondname,status='old')
 ! Get parameter filenames and open the parameter files.
       call getenv_loc('BONDPAR',bondname)
       open (ibond,file=bondname,status='old')

+      call getenv_loc('BONDPAR_NUCL',bondname_nucl)
+      open (ibond_nucl,file=bondname_nucl,status='old')
+

       call getenv_loc('THETPAR',thetname)
       open (ithep,file=thetname,status='old')
       call getenv_loc('ROTPAR',rotname)
       call getenv_loc('THETPAR',thetname)
       open (ithep,file=thetname,status='old')
       call getenv_loc('ROTPAR',rotname)

@@ -3817,10 +4667,23 @@
       open (ielep,file=elename,status='old')
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old')
       open (ielep,file=elename,status='old')
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old')

+
+      open (ithep_nucl,file=thetname_nucl,status='old')
+      call getenv_loc('ROTPAR_NUCL',rotname_nucl)
+      open (irotam_nucl,file=rotname_nucl,status='old')
+      call getenv_loc('TORPAR_NUCL',torname_nucl)
+      open (itorp_nucl,file=torname_nucl,status='old')
+      call getenv_loc('TORDPAR_NUCL',tordname_nucl)
+      open (itordp_nucl,file=tordname_nucl,status='old')
+      call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
+      open (isidep_nucl,file=sidename_nucl,status='old')
+

       call getenv_loc('LIPTRANPAR',liptranname)
       open (iliptranpar,file=liptranname,status='old')
       call getenv_loc('TUBEPAR',tubename)
       open (itube,file=tubename,status='old')
       call getenv_loc('LIPTRANPAR',liptranname)
       open (iliptranpar,file=liptranname,status='old')
       call getenv_loc('TUBEPAR',tubename)
       open (itube,file=tubename,status='old')

+      call getenv_loc('IONPAR',ionname)
+      open (iion,file=ionname,status='old')

 #else
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
         readonly)
 #else
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
         readonly)

@@ -3829,6 +4692,8 @@
 ! Get parameter filenames and open the parameter files.
       call getenv_loc('BONDPAR',bondname)
       open (ibond,file=bondname,status='old',action='read')
 ! Get parameter filenames and open the parameter files.
       call getenv_loc('BONDPAR',bondname)
       open (ibond,file=bondname,status='old',action='read')

+      call getenv_loc('BONDPAR_NUCL',bondname_nucl)
+      open (ibond_nucl,file=bondname_nucl,status='old',action='read')

       call getenv_loc('THETPAR',thetname)
       open (ithep,file=thetname,status='old',action='read')
       call getenv_loc('ROTPAR',rotname)
       call getenv_loc('THETPAR',thetname)
       open (ithep,file=thetname,status='old',action='read')
       call getenv_loc('ROTPAR',rotname)

@@ -3851,16 +4716,50 @@
       open (ielep,file=elename,status='old',readonly)
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old',readonly)
       open (ielep,file=elename,status='old',readonly)
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old',readonly)

+
+      call getenv_loc('THETPAR_NUCL',thetname_nucl)
+      open (ithep_nucl,file=thetname_nucl,status='old',action='read')
+      call getenv_loc('ROTPAR_NUCL',rotname_nucl)
+      open (irotam_nucl,file=rotname_nucl,status='old',action='read')
+      call getenv_loc('TORPAR_NUCL',torname_nucl)
+      open (itorp_nucl,file=torname_nucl,status='old',action='read')
+      call getenv_loc('TORDPAR_NUCL',tordname_nucl)
+      open (itordp_nucl,file=tordname_nucl,status='old',action='read')
+      call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
+      open (isidep_nucl,file=sidename_nucl,status='old',action='read')
+      call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
+      open (isidep_scbase,file=sidename_scbase,status='old',action='read')
+      call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
+      open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
+      call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
+      open (isidep_scpho,file=pepname_scpho,status='old',action='read')
+      call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
+      open (isidep_peppho,file=pepname_peppho,status='old',action='read')
+
+

       call getenv_loc('LIPTRANPAR',liptranname)
       open (iliptranpar,file=liptranname,status='old',action='read')
       call getenv_loc('TUBEPAR',tubename)
       open (itube,file=tubename,status='old',action='read')
       call getenv_loc('LIPTRANPAR',liptranname)
       open (iliptranpar,file=liptranname,status='old',action='read')
       call getenv_loc('TUBEPAR',tubename)
       open (itube,file=tubename,status='old',action='read')

+      call getenv_loc('IONPAR',ionname)
+      open (iion,file=ionname,status='old',action='read')

 
 #ifndef CRYST_SC
       call getenv_loc('ROTPARPDB',rotname_pdb)
       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
 #endif
 #endif
 
 #ifndef CRYST_SC
       call getenv_loc('ROTPARPDB',rotname_pdb)
       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
 #endif
 #endif

+      call getenv_loc('SCPPAR_NUCL',scpname_nucl)
+#if defined(WINIFL) || defined(WINPGI)
+      open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
+#elif (defined CRAY)  || (defined AIX)
+      open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
+#elif (defined G77)
+      open (iscpp_nucl,file=scpname_nucl,status='old')
+#else
+      open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
+#endif
+

 #ifndef OLDSCP
 !
 ! 8/9/01 In the newest version SCp interaction constants are read from a file
 #ifndef OLDSCP
 !
 ! 8/9/01 In the newest version SCp interaction constants are read from a file

@@ -4052,6 +4951,7 @@
 
       open(irest2,file=rest2name,status='unknown')
       read(irest2,*) totT,EK,potE,totE,t_bath
 
       open(irest2,file=rest2name,status='unknown')
       read(irest2,*) totT,EK,potE,totE,t_bath

+      totTafm=totT

 !      do i=1,2*nres
 ! AL 4/17/17: Now reading d_t(0,:) too
       do i=0,2*nres
 !      do i=1,2*nres
 ! AL 4/17/17: Now reading d_t(0,:) too
       do i=0,2*nres