Merge branch 'UCGM' of mmka.chem.univ.gda.pl:unres4 into UCGM
[unres4.git] / source / unres / io_config.f90
index e6b7f0a..87ddf94 100644 (file)
 ! 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) ='-'
       
       ii=0
       do i=1,nres
+         ftors(i)=ftors(1)
         if ( secstruc(i) .eq. 'H') then
 ! Helix restraints for this residue               
            ii=ii+1 
       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
-      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(nbondterm(ntyp)) !(ntyp)
       allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
       allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
-      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))
       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
         endif
       enddo
 #else
+      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),&
           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)
          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.
 
 #endif      
       if (lprint) then
+        l=3
         write (iout,'(/a/)') 'Torsional constants:'
         do i=1,nsccortyp
           do j=1,nsccortyp
 
       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
-          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
-        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
-        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
-        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
-        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
-        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
 !---------------------- 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
          write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
                              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)
 ! 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(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),&
       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
+        endif
       acavtub(:)=0.0d0
       bcavtub(:)=0.0d0
       ccavtub(:)=0.0d0
           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)
           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   ',&
          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
          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)
           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)
        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
       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
       use control_data
       use compare_data
       use MPI_data
-      use control, only: rescode,sugarcode
+!      use control, only: rescode,sugarcode
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.LOCAL'
       character(len=80) :: card
       real(kind=8),dimension(3,20) :: sccor
       integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin      !rescode,
-      integer :: isugar,molecprev
+      integer :: isugar,molecprev,firstion
       character*1 :: sugar
       real(kind=8) :: cou
       real(kind=8),dimension(3) :: ccc
 !      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
+      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)
-
+      if(.not. allocated(istype)) allocate(istype(maxres))
       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 (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)
+              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)
             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
 !              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. &
               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 (molecule.ne.5) molecprev=molecule
            molecule=5
            nres_molec(molecule)=nres_molec(molecule)+1
-           itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
+           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
 ! Calculate dummy residue coordinates inside the "chain" of a multichain
 ! system
       nres=ires
-      if (ires_old.ne.ires) nres_molec(molecule)=nres_molec(molecule)-2
-!      print *,'I have', nres_molec(:)
+      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
 ! 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
-           if (unres_pdb) then
+!           if (unres_pdb) then
 ! 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
+!            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
+!             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
-          endif   !unres_pdb
+!          endif   !unres_pdb
          else     !itype(i+1,1).eq.ntyp1
-          if (unres_pdb) then
+!          if (unres_pdb) then
 ! 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
-          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
-          endif !unres_pdb
+!          endif !unres_pdb
          endif !itype(i+1,1).eq.ntyp1
         endif  !itype.eq.ntyp1
 
           if (molecule.eq.5) molecule=molecprev
         itype(nres,molecule)=ntyp1_molec(molecule)
         nres_molec(molecule)=nres_molec(molecule)+1
-        if (unres_pdb) then
+!        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)
-          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
-        endif
+!        endif
       endif
 !     print *,'I have',nres, nres_molec(:)
 
       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
-        write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.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
        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)') &
+        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
 ! NOW LETS ROCK! SORTING
       allocate(c_temporary(3,2*nres))
       allocate(itype_temporary(nres,5))
-      allocate(molnum(nres))
+      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
 
           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
          istype(i)=istype_temp(i)
         enddo
        enddo
-      if (itype(1,1).eq.ntyp1) then
-        nsup=nsup-1
-        nstart_sup=2
-        if (unres_pdb) then
+!      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
+!          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))') &
+      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),i3,3f8.3,5x,3f8.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
         write (iout,*) "symetr", symetr
       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 ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
       chain_length=lll-1
 !       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
         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
             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(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')
-  110 format (a,'(',i3,')',6f12.5)
+  110 format (a,'(',i5,')',6f12.5)
      
       enddo
 !c enddiag
       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,'(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      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
-      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 
-      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
 !      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.)
       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(me.eq.king.or..not.out1file) &
          write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
           eta
-        gamp=scal_fric*(pstok(i)+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,1)+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,1),&
-           gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
+           gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
          enddo
         endif
       else if (tbf) then
       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" 
       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)
       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('IONPAR',ionname)
+      open (iion,file=ionname,status='old',action='read')
 
 #ifndef CRYST_SC
       call getenv_loc('ROTPARPDB',rotname_pdb)