working side-chain phosphate
[unres4.git] / source / unres / io_config.f90
index 277b6ba..2d4a625 100644 (file)
       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,&
 !
       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(msc(ntyp+1)) !(ntyp+1)
-      allocate(isc(ntyp+1)) !(ntyp+1)
-      allocate(restok(ntyp+1)) !(ntyp+1)
+      allocate(msc(ntyp+1,5)) !(ntyp+1)
+      allocate(isc(ntyp+1,5)) !(ntyp+1)
+      allocate(restok(ntyp+1,5)) !(ntyp+1)
       allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
       allocate(long_r_sidechain(ntyp))
       allocate(short_r_sidechain(ntyp))
         endif
       enddo
 #else
-      read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp,ip,pstok
-      do i=1,ntyp
+      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),&
-         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
         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'
-        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
           write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
-            vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(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 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)
       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)
          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.
       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
       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))
       end select
       continue
       close (isidep)
+
 !-----------------------------------------------------------------------
 ! Calculate the "working" parameters of SC interactions.
 
          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
        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(j,i)
+        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)))
+
+!      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)
+       
+       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
+
 
       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
       bad(:,:)=0.0D0
       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
 !
 ! 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, only: rescode
+      use control, only: rescode,sugarcode
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.LOCAL'
 !      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
       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),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)
-
+      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
+      molecule=1
+      counter=0
 !      write (2,*) "UNRES_PDB",unres_pdb
       ires=0
       ires_old=0
         else if (card(:3).eq.'TER') then
 ! End current chain
           ires_old=ires+2
+          ishift=ishift+1
           ishift1=ishift1+1
-          itype(ires_old,1)=ntyp1
-          itype(ires_old-1,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
 ! 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 (res.ne.'GLY' .and. res.ne. 'ACE') then
                 ishift=ishift-1
                 itype(1,1)=ntyp1
+                nres_molec(molecule)=nres_molec(molecule)+1
               endif
               ires=ires-ishift+ishift1
               ires_old=ires
               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
-            endif
+            endif 
+!            print *,'atom',ires,atom
             if (res.eq.'ACE' .or. res.eq.'NHE') then
               itype(ires,1)=10
             else
-              itype(ires,1)=rescode(ires,res,0,1)
+             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
+             molecule=2
+              itype(ires,molecule)=rescode(ires,res(2:4),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
           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)') &
             ires,itype(ires,1),res,(c(j,ires),j=1,3)
 #endif
             iii=iii+1
+              nres_molec(molecule)=nres_molec(molecule)+1
             do j=1,3
               sccor(j,iii)=c(j,ires)
             enddo
+          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
+              read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+            else
+              iii=iii+1
+              read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+            endif
 !            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. &
-                   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,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)
+              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)
+           itype(ires,molecule)=rescode(ires,res(2:4),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
+      if ((ires_old.ne.ires).and.(molecule.ne.5)) &
+         nres_molec(molecule)=nres_molec(molecule)-2
+!      print *,'I have', nres_molec(:)
+      
+      do k=1,4 ! ions are without dummy 
+       if (nres_molec(k).eq.0) cycle
       do i=2,nres-1
 !        write (iout,*) i,itype(i,1)
 !        if (itype(i,1).eq.ntyp1) then
 !            dc(j,i)=c(j,i)
 !          enddo
 !        endif
-        if (itype(i,1).eq.ntyp1) then
-         if (itype(i+1,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
 ! 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
               e2(2)=1.0d0
               e2(3)=0.0d0
             endif !fail
-            print *,i,'a tu?'
+!            print *,i,'a tu?'
             do j=1,3
              c(j,i)=c(j,i-1)-1.9d0*e2(j)
             enddo
         endif  !itype.eq.ntyp1
 
       enddo
+      enddo
 ! Calculate the CM of the last side chain.
       if (iii.gt.0)  then
       if (unres_pdb) then
 !      nres=ires
       nsup=nres
       nstart_sup=1
-      if (itype(nres,1).ne.10) then
+!      print *,"molecule",molecule
+      if ((itype(nres,1).ne.10)) then
         nres=nres+1
-        itype(nres,1)=ntyp1
+          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
           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
         enddo
         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
         c(j,nres+1)=c(j,1)
         c(j,2*nres)=c(j,nres)
       enddo
+      
       if (itype(1,1).eq.ntyp1) then
         nsup=nsup-1
         nstart_sup=2
         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
       write (iout,'(a,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,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
-          restyp(itype(ires,1),1),ires,(c(j,ires),j=1,3),&
+        write (iout,'(5(a3,1x),i3,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,nres+ires),j=1,3)
        enddo
       endif
+! NOW LETS ROCK! SORTING
+      allocate(c_temporary(3,2*nres))
+      allocate(itype_temporary(nres,5))
+      allocate(molnum(nres))
+      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,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),i3,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
       kkk=1
       lll=0
       cou=1
+        write (iout,*) "symetr", symetr
       do i=1,nres
       lll=lll+1
 !c      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
 !      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.)
         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
-          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
+
         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
-          write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),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
 ! 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)
       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')
 ! 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 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('IONPAR',ionname)
+      open (iion,file=ionname,status='old',action='read')
 
 !      print *,"Processor",myrank," opened file ISIDEP" 
 !      print *,"Processor",myrank," opened parameter files" 
 ! 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)
       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('IONPAR',ionname)
+      open (iion,file=ionname,status='old')
 #else
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
         readonly)
 ! 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)
       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('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
+      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