distance constraint energy calculation in wham
[unres4.git] / source / wham / io_wham.F90
index 958e244..187502e 100644 (file)
@@ -51,6 +51,8 @@
 ! Get parameter filenames and open the parameter files.
       call mygetenv('BONDPAR',bondname)
       open (ibond,file=bondname,status='old')
+      call mygetenv('BONDPAR_NUCL',bondname_nucl)
+      open (ibond_nucl,file=bondname_nucl,status='old')
       call mygetenv('THETPAR',thetname)
       open (ithep,file=thetname,status='old')
       call mygetenv('ROTPAR',rotname)
       open (isidep,file=sidename,status='old')
       call mygetenv('SIDEP',sidepname)
       open (isidep1,file=sidepname,status="old")
+      call mygetenv('THETPAR_NUCL',thetname_nucl)
+      open (ithep_nucl,file=thetname_nucl,status='old')
+      call mygetenv('ROTPAR_NUCL',rotname_nucl)
+      open (irotam_nucl,file=rotname_nucl,status='old')
+      call mygetenv('TORPAR_NUCL',torname_nucl)
+      open (itorp_nucl,file=torname_nucl,status='old')
+      call mygetenv('TORDPAR_NUCL',tordname_nucl)
+      open (itordp_nucl,file=tordname_nucl,status='old')
+      call mygetenv('SIDEPAR_NUCL',sidename_nucl)
+      open (isidep_nucl,file=sidename_nucl,status='old')
+      call mygetenv('SIDEPAR_SCBASE',sidename_scbase)
+      open (isidep_scbase,file=sidename_scbase,status='old')
+      call mygetenv('PEPPAR_PEPBASE',pepname_pepbase)
+      open (isidep_pepbase,file=pepname_pepbase,status='old')
+      call mygetenv('SCPAR_PHOSPH',pepname_scpho)
+      open (isidep_scpho,file=pepname_scpho,status='old')
+      call mygetenv('PEPPAR_PHOSPH',pepname_peppho)
+      open (isidep_peppho,file=pepname_peppho,status='old')
+      call mygetenv('LIPTRANPAR',liptranname)
+      open (iliptranpar,file=liptranname,status='old',action='read')
+      call mygetenv('TUBEPAR',tubename)
+      open (itube,file=tubename,status='old',action='read')
+      call mygetenv('IONPAR',ionname)
+      open (iion,file=ionname,status='old',action='read')
+
+      call mygetenv('SCPPAR_NUCL',scpname_nucl)
+      open (iscpp_nucl,file=scpname_nucl,status='old')
+
+
 #ifndef OLDSCP
 !
 ! 8/9/01 In the newest version SCp interaction constants are read from a file
       integer :: iparm,nkcctyp
 !el      real(kind=8) :: ip,mp
       real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,epsijlip,&
-                sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm
+                sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,    &
+                epspeptube,epssctube,sigmapeptube,      &
+                sigmasctube
+
       real(kind=8) :: v0ij,v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,rjunk,&
                 res1
       integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n,jj
@@ -427,11 +461,39 @@ allocate(ww(max_eneW))
       wscloc=ww(12)
       wtor=ww(13)
       wtor_d=ww(14)
+      wstrain=ww(15)
       wvdwpp=ww(16)
       wbond=ww(18)
       wsccor=ww(19)
       wcatcat=ww(42)
       wcatprot=ww(41)
+      wcorr3_nucl=ww(38)
+      wcorr_nucl=ww(37)
+      wtor_d_nucl=ww(36)
+      wtor_nucl=ww(35)
+      wsbloc=ww(34)
+      wang_nucl=ww(33)
+      wbond_nucl=ww(32)
+      welsb=ww(31)
+      wvdwsb=ww(30)
+      welpsb=ww(29)
+      wvdwpsb=ww(28)
+      welpp=ww(27)
+      wvdwpp_nucl=ww(26)
+      wscbase=ww(46)
+      wpepbase=ww(47)
+      wscpho=ww(48)
+      wpeppho=ww(49)
+!      print *,"KURWA",ww(48)
+!        "WSCBASE   ","WPEPBASE  ","WSCPHO    ","WPEPPHO   "
+!        "WVDWPP    ","WELPP     ","WVDWPSB   ","WELPSB    ","WVDWSB    ",&
+!        "WELSB     ","WBOND     ","WANG      ","WSBLOC    ","WTOR      ",&
+!        "WTORD     ","WCORR     ","WCORR3    ","WNULL     ","WNULL     ",&
+!        "WCATPROT  ","WCATCAT  
+!       +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+!       +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+!       +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+!       +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
 
       endif
 !
@@ -451,14 +513,32 @@ allocate(ww(max_eneW))
       weights(12)=wscloc
       weights(13)=wtor
       weights(14)=wtor_d
-      weights(15)=0 !wstrain !
-      weights(16)=0 !wvdwpp !
+      weights(15)=wstrain !0
+      weights(16)=wvdwpp !
       weights(17)=wbond
       weights(18)=0 !scal14 !
       weights(21)=wsccor
       weights(42)=wcatprot
       weights(41)=wcatcat
-
+      weights(26)=    wvdwpp_nucl 
+
+      weights(27) =welpp  
+      weights(28) =wvdwpsb
+      weights(29) =welpsb 
+      weights(30) =wvdwsb 
+      weights(31) =welsb  
+      weights(32) =wbond_nucl  
+      weights(33) =wang_nucl   
+      weights(34) =wsbloc 
+      weights(35) =wtor_nucl   
+      weights(36) =wtor_d_nucl 
+      weights(37) =wcorr_nucl  
+      weights(38) =wcorr3_nucl 
+      weights(41) =wcatcat
+      weights(42) =wcatprot
+      weights(46) =wscbase
+      weights(47) =wscpho
+      weights(48) =wpeppho
 ! el--------
       call card_concat(controlcard,.false.)
 
@@ -536,6 +616,10 @@ allocate(ww(max_eneW))
       allocate(nbondterm(ntyp)) !(ntyp)
       allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
       allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
+      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)
+
 !el      allocate(msc(ntyp+1)) !(ntyp+1)
 !el      allocate(isc(ntyp+1)) !(ntyp+1)
 !el      allocate(restok(ntyp+1)) !(ntyp+1)
@@ -595,6 +679,17 @@ allocate(ww(max_eneW))
             read (iion,*)  (catprm(i,k),i=1,ncatprotparm)
             enddo
             print *, catprm
+      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),rjunk,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
 
 
 !----------------------------------------------------
@@ -933,6 +1028,67 @@ allocate(ww(max_eneW))
       ENDIF ! TOR_MODE
 
 #endif
+!--------------- Reading theta parameters for nucleic acid-------
+      read (ithep_nucl,*) 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,*) (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)') t1,t2,t3
+            read (ithep_nucl,*) aa0thet_nucl(i,j,k)
+            read (ithep_nucl,*)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
+            read (ithep_nucl,*) &
+            (((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,*) &
+           (((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(bsc(maxlob,ntyp)) !(maxlob,ntyp)
@@ -1070,6 +1226,24 @@ allocate(ww(max_eneW))
       enddo
 #endif
       close(irotam)
+!---------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,*)
+        do j=1,9
+          read(irotam_nucl,*) 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 (ifourier,*) nloctyp
 !el write(iout,*)"nloctyp",nloctyp
@@ -1651,6 +1825,11 @@ allocate(ww(max_eneW))
       enddo
       enddo
       endif
+#ifndef NEWCORR
+      do i=1,ntyp1
+        itype2loc(i)=itortyp(i)
+      enddo
+#endif
       ELSE IF (TOR_MODE.eq.1) THEN
 
 !C read valence-torsional parameters
@@ -2117,6 +2296,10 @@ allocate(ww(max_eneW))
       allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
       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),&
+        tubetranene(ntyp1))
       do i=1,ntyp1
         do j=1,ntyp1
           aa_aq(i,j)=0.0D0
@@ -2215,6 +2398,271 @@ allocate(ww(max_eneW))
          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
+      sigmapeptube=sigmapeptube**6
+      sigeps=dsign(1.0D0,epspeptube)
+      epspeptube=dabs(epspeptube)
+      pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
+      pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
+      write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
+      do i=1,ntyp
+       read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
+      ccavtub(i),dcavtub(i),tubetranene(i)
+       sigmasctube=sigmasctube**6
+       sigeps=dsign(1.0D0,epssctube)
+       epssctube=dabs(epssctube)
+       sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
+       sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
+      write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
+      enddo
+!-----------------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------------------------------------- 
+!--------------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)
       do i=1,ntyp
@@ -2233,6 +2681,48 @@ allocate(ww(max_eneW))
         enddo
       enddo
 #endif
+      allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
+      read (itorp_nucl,*) 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,*) &
+        (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,*) 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,*) 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_nucl,*) 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
+
 
 !elwrite(iout,*) "parmread kontrol before oldscp"
 !
@@ -2283,6 +2773,20 @@ allocate(ww(max_eneW))
         enddo
       endif
 #endif
+      allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
+
+      do i=1,ntyp_molec(2)
+        read (iscpp_nucl,*) 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
 !
@@ -3109,6 +3613,7 @@ allocate(ww(max_eneW))
                 do k=1,3
                   cref(k,j+i,kkk)=cref(k,j,kkk)
                 enddo
+                write(iout,*) "J",j,"J+I",j+i
                 phi_ref(j+i)=phi_ref(j)
                 theta_ref(j+i)=theta_ref(j)
                 alph_ref(j+i)=alph_ref(j)
@@ -3223,8 +3728,9 @@ allocate(ww(max_eneW))
 !      include 'COMMON.HEADER'
 !      include 'COMMON.SBRIDGE'
       character(len=50) :: tytul
-      character(len=1),dimension(10) :: chainid=reshape((/'A','B','C',&
-                      'D','E','F','G','H','I','J'/),shape(chainid))
+      character(len=1),dimension(24) :: chainid=reshape((/'A','B','C',&
+                      'D','E','F','G','H','I','J','K','L','M','N','O',&
+                      'P','Q','R','S','V','W','X','Y','Z'/),shape(chainid))
       integer,dimension(nres) :: ica !(maxres)
       real(kind=8) :: temp,efree,etot,entropy,rmsdev
       integer :: ii,i,j,iti,ires,iatom,ichain,mnum
@@ -3241,9 +3747,11 @@ allocate(ww(max_eneW))
         mnum=molnum(i)
         iti=itype(i,mnum)
         if (iti.eq.ntyp1) then
+          if (itype(i-1,molnum(i-1)).eq.ntyp1) then
           ichain=ichain+1
           ires=0
           write (ipdb,'(a)') 'TER'
+          endif
         else
         ires=ires+1
         iatom=iatom+1