projects / unres4.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
Merge branch 'UCGM' of mmka.chem.univ.gda.pl:unres4 into UCGM
[unres4.git] / source / unres / io_config.f90
diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90
index 6186282..87ddf94 100644 (file)
--- a/source/unres/io_config.f90
+++ b/source/unres/io_config.f90
@@ -529,9 +529,9 @@
 ! read secondary structure prediction from JPRED here!
 !      read(isecpred,'(A80)',err=100,end=100) line
 !      read(line,'(f10.3)',err=110) ftors
 ! read secondary structure prediction from JPRED here!
 !      read(isecpred,'(A80)',err=100,end=100) line
 !      read(line,'(f10.3)',err=110) ftors
-       read(isecpred,'(f10.3)',err=110) ftors
+       read(isecpred,'(f10.3)',err=110) ftors(1)
 
 
-      write (iout,*) 'FTORS factor =',ftors
+      write (iout,*) 'FTORS factor =',ftors(1)
 ! initialize secstruc to any
        do i=1,nres
         secstruc(i) ='-'
 ! initialize secstruc to any
        do i=1,nres
         secstruc(i) ='-'
@@ -553,6 +553,7 @@
       
       ii=0
       do i=1,nres
       
       ii=0
       do i=1,nres
+         ftors(i)=ftors(1)
         if ( secstruc(i) .eq. 'H') then
 ! Helix restraints for this residue               
            ii=ii+1 
         if ( secstruc(i) .eq. 'H') then
 ! Helix restraints for this residue               
            ii=ii+1 
@@ -735,7 +736,7 @@
       character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
       logical :: lprint,LaTeX
       real(kind=8),dimension(3,3,maxlob) :: blower     !(3,3,maxlob)
       character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
       logical :: lprint,LaTeX
       real(kind=8),dimension(3,3,maxlob) :: blower     !(3,3,maxlob)
-      real(kind=8),dimension(13) :: b
+      real(kind=8),dimension(13) :: buse
       character(len=3) :: lancuch      !,ucase
 !el  local variables
       integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
       character(len=3) :: lancuch      !,ucase
 !el  local variables
       integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
@@ -782,9 +783,6 @@
       allocate(nbondterm(ntyp)) !(ntyp)
       allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
       allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
       allocate(nbondterm(ntyp)) !(ntyp)
       allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
       allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
-      allocate(msc(ntyp+1,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))
       allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
       allocate(long_r_sidechain(ntyp))
       allocate(short_r_sidechain(ntyp))
@@ -792,6 +790,10 @@
       dsc_inv(:)=0.0d0
 
 #ifdef CRYST_BOND
       dsc_inv(:)=0.0d0
 
 #ifdef CRYST_BOND
+      allocate(msc(ntyp+1)) !(ntyp+1)
+      allocate(isc(ntyp+1)) !(ntyp+1)
+      allocate(restok(ntyp+1)) !(ntyp+1)
+
       read (ibond,*) vbldp0,akp,mp,ip,pstok
       do i=1,ntyp
         nbondterm(i)=1
       read (ibond,*) vbldp0,akp,mp,ip,pstok
       do i=1,ntyp
         nbondterm(i)=1
@@ -804,6 +806,15 @@
         endif
       enddo
 #else
         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),&
       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),&
@@ -1897,6 +1908,7 @@
 
 #endif      
       if (lprint) then
 
 #endif      
       if (lprint) then
+        l=3
         write (iout,'(/a/)') 'Torsional constants:'
         do i=1,nsccortyp
           do j=1,nsccortyp
         write (iout,'(/a/)') 'Torsional constants:'
         do i=1,nsccortyp
           do j=1,nsccortyp
@@ -1941,88 +1953,88 @@
 
       do i=0,nloctyp-1
         read (ifourier,*,end=115,err=115)
 
       do i=0,nloctyp-1
         read (ifourier,*,end=115,err=115)
-        read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
+        read (ifourier,*,end=115,err=115) (buse(ii),ii=1,13)
         if (lprint) then
           write (iout,*) 'Type',i
         if (lprint) then
           write (iout,*) 'Type',i
-          write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
+          write (iout,'(a,i2,a,f10.5)') ('buse(',ii,')=',buse(ii),ii=1,13)
         endif
         endif
-        B1(1,i)  = b(3)
-        B1(2,i)  = b(5)
-        B1(1,-i) = b(3)
-        B1(2,-i) = -b(5)
-!        b1(1,i)=0.0d0
-!        b1(2,i)=0.0d0
-        B1tilde(1,i) = b(3)
-        B1tilde(2,i) =-b(5)
-        B1tilde(1,-i) =-b(3)
-        B1tilde(2,-i) =b(5)
-!        b1tilde(1,i)=0.0d0
-!        b1tilde(2,i)=0.0d0
-        B2(1,i)  = b(2)
-        B2(2,i)  = b(4)
-        B2(1,-i)  =b(2)
-        B2(2,-i)  =-b(4)
-
-!        b2(1,i)=0.0d0
-!        b2(2,i)=0.0d0
-        CC(1,1,i)= b(7)
-        CC(2,2,i)=-b(7)
-        CC(2,1,i)= b(9)
-        CC(1,2,i)= b(9)
-        CC(1,1,-i)= b(7)
-        CC(2,2,-i)=-b(7)
-        CC(2,1,-i)=-b(9)
-        CC(1,2,-i)=-b(9)
+        B1(1,i)  = buse(3)
+        B1(2,i)  = buse(5)
+        B1(1,-i) = buse(3)
+        B1(2,-i) = -buse(5)
+!        buse1(1,i)=0.0d0
+!        buse1(2,i)=0.0d0
+        B1tilde(1,i) = buse(3)
+        B1tilde(2,i) =-buse(5)
+        B1tilde(1,-i) =-buse(3)
+        B1tilde(2,-i) =buse(5)
+!        buse1tilde(1,i)=0.0d0
+!        buse1tilde(2,i)=0.0d0
+        B2(1,i)  = buse(2)
+        B2(2,i)  = buse(4)
+        B2(1,-i)  =buse(2)
+        B2(2,-i)  =-buse(4)
+
+!        buse2(1,i)=0.0d0
+!        buse2(2,i)=0.0d0
+        CC(1,1,i)= buse(7)
+        CC(2,2,i)=-buse(7)
+        CC(2,1,i)= buse(9)
+        CC(1,2,i)= buse(9)
+        CC(1,1,-i)= buse(7)
+        CC(2,2,-i)=-buse(7)
+        CC(2,1,-i)=-buse(9)
+        CC(1,2,-i)=-buse(9)
 !        CC(1,1,i)=0.0d0
 !        CC(2,2,i)=0.0d0
 !        CC(2,1,i)=0.0d0
 !        CC(1,2,i)=0.0d0
 !        CC(1,1,i)=0.0d0
 !        CC(2,2,i)=0.0d0
 !        CC(2,1,i)=0.0d0
 !        CC(1,2,i)=0.0d0
-        Ctilde(1,1,i)=b(7)
-        Ctilde(1,2,i)=b(9)
-        Ctilde(2,1,i)=-b(9)
-        Ctilde(2,2,i)=b(7)
-        Ctilde(1,1,-i)=b(7)
-        Ctilde(1,2,-i)=-b(9)
-        Ctilde(2,1,-i)=b(9)
-        Ctilde(2,2,-i)=b(7)
+        Ctilde(1,1,i)=buse(7)
+        Ctilde(1,2,i)=buse(9)
+        Ctilde(2,1,i)=-buse(9)
+        Ctilde(2,2,i)=buse(7)
+        Ctilde(1,1,-i)=buse(7)
+        Ctilde(1,2,-i)=-buse(9)
+        Ctilde(2,1,-i)=buse(9)
+        Ctilde(2,2,-i)=buse(7)
 
 !        Ctilde(1,1,i)=0.0d0
 !        Ctilde(1,2,i)=0.0d0
 !        Ctilde(2,1,i)=0.0d0
 !        Ctilde(2,2,i)=0.0d0
 
 !        Ctilde(1,1,i)=0.0d0
 !        Ctilde(1,2,i)=0.0d0
 !        Ctilde(2,1,i)=0.0d0
 !        Ctilde(2,2,i)=0.0d0
-        DD(1,1,i)= b(6)
-        DD(2,2,i)=-b(6)
-        DD(2,1,i)= b(8)
-        DD(1,2,i)= b(8)
-        DD(1,1,-i)= b(6)
-        DD(2,2,-i)=-b(6)
-        DD(2,1,-i)=-b(8)
-        DD(1,2,-i)=-b(8)
+        DD(1,1,i)= buse(6)
+        DD(2,2,i)=-buse(6)
+        DD(2,1,i)= buse(8)
+        DD(1,2,i)= buse(8)
+        DD(1,1,-i)= buse(6)
+        DD(2,2,-i)=-buse(6)
+        DD(2,1,-i)=-buse(8)
+        DD(1,2,-i)=-buse(8)
 !        DD(1,1,i)=0.0d0
 !        DD(2,2,i)=0.0d0
 !        DD(2,1,i)=0.0d0
 !        DD(1,2,i)=0.0d0
 !        DD(1,1,i)=0.0d0
 !        DD(2,2,i)=0.0d0
 !        DD(2,1,i)=0.0d0
 !        DD(1,2,i)=0.0d0
-        Dtilde(1,1,i)=b(6)
-        Dtilde(1,2,i)=b(8)
-        Dtilde(2,1,i)=-b(8)
-        Dtilde(2,2,i)=b(6)
-        Dtilde(1,1,-i)=b(6)
-        Dtilde(1,2,-i)=-b(8)
-        Dtilde(2,1,-i)=b(8)
-        Dtilde(2,2,-i)=b(6)
+        Dtilde(1,1,i)=buse(6)
+        Dtilde(1,2,i)=buse(8)
+        Dtilde(2,1,i)=-buse(8)
+        Dtilde(2,2,i)=buse(6)
+        Dtilde(1,1,-i)=buse(6)
+        Dtilde(1,2,-i)=-buse(8)
+        Dtilde(2,1,-i)=buse(8)
+        Dtilde(2,2,-i)=buse(6)
 
 !        Dtilde(1,1,i)=0.0d0
 !        Dtilde(1,2,i)=0.0d0
 !        Dtilde(2,1,i)=0.0d0
 !        Dtilde(2,2,i)=0.0d0
 
 !        Dtilde(1,1,i)=0.0d0
 !        Dtilde(1,2,i)=0.0d0
 !        Dtilde(2,1,i)=0.0d0
 !        Dtilde(2,2,i)=0.0d0
-        EE(1,1,i)= b(10)+b(11)
-        EE(2,2,i)=-b(10)+b(11)
-        EE(2,1,i)= b(12)-b(13)
-        EE(1,2,i)= b(12)+b(13)
-        EE(1,1,-i)= b(10)+b(11)
-        EE(2,2,-i)=-b(10)+b(11)
-        EE(2,1,-i)=-b(12)+b(13)
-        EE(1,2,-i)=-b(12)-b(13)
+        EE(1,1,i)= buse(10)+buse(11)
+        EE(2,2,i)=-buse(10)+buse(11)
+        EE(2,1,i)= buse(12)-buse(13)
+        EE(1,2,i)= buse(12)+buse(13)
+        EE(1,1,-i)= buse(10)+buse(11)
+        EE(2,2,-i)=-buse(10)+buse(11)
+        EE(2,1,-i)=-buse(12)+buse(13)
+        EE(1,2,-i)=-buse(12)-buse(13)
 
 !        ee(1,1,i)=1.0d0
 !        ee(2,2,i)=1.0d0
 
 !        ee(1,1,i)=1.0d0
 !        ee(2,2,i)=1.0d0
@@ -2146,6 +2158,8 @@
 !---------------------- GB or BP potential -----------------------------
        case(3:4)
 !   30 do i=1,ntyp
 !---------------------- GB or BP potential -----------------------------
        case(3:4)
 !   30 do i=1,ntyp
+!        print *,"I AM in SCELE",scelemode
+        if (scelemode.eq.0) then
         do i=1,ntyp
          read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
         enddo
         do i=1,ntyp
          read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
         enddo
@@ -2173,6 +2187,81 @@
          write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
                              chip(i),alp(i),i=1,ntyp)
         endif
          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)
 !      goto 50
 !--------------------- GBV potential -----------------------------------
        case(5)
@@ -2201,9 +2290,11 @@
 ! Calculate the "working" parameters of SC interactions.
 
 !el from module energy - COMMON.INTERACT-------
 ! Calculate the "working" parameters of SC interactions.
 
 !el from module energy - COMMON.INTERACT-------
-      allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
+      allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
+      if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
       allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
       allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
-      allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
+      if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
+      allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
       allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
         dcavtub(ntyp1))
       allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
       allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
         dcavtub(ntyp1))
       allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
@@ -2212,9 +2303,11 @@
       bb_aq(:,:)=0.0D0
       aa_lip(:,:)=0.0D0
       bb_lip(:,:)=0.0D0
       bb_aq(:,:)=0.0D0
       aa_lip(:,:)=0.0D0
       bb_lip(:,:)=0.0D0
+         if (scelemode.eq.0) then
       chi(:,:)=0.0D0
       sigma(:,:)=0.0D0
       r0(:,:)=0.0D0
       chi(:,:)=0.0D0
       sigma(:,:)=0.0D0
       r0(:,:)=0.0D0
+        endif
       acavtub(:)=0.0d0
       bcavtub(:)=0.0d0
       ccavtub(:)=0.0d0
       acavtub(:)=0.0d0
       bcavtub(:)=0.0d0
       ccavtub(:)=0.0d0
@@ -2230,6 +2323,7 @@
           epslip(i,j)=epslip(j,i)
         enddo
       enddo
           epslip(i,j)=epslip(j,i)
         enddo
       enddo
+         if (scelemode.eq.0) then
       do i=1,ntyp
         do j=i,ntyp
           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
       do i=1,ntyp
         do j=i,ntyp
           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
@@ -2238,6 +2332,7 @@
           rs0(j,i)=rs0(i,j)
         enddo
       enddo
           rs0(j,i)=rs0(i,j)
         enddo
       enddo
+      endif
       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
        'Working parameters of the SC interactions:',&
        '     a    ','     b    ','   augm   ','  sigma ','   r0   ',&
       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
        'Working parameters of the SC interactions:',&
        '     a    ','     b    ','   augm   ','  sigma ','   r0   ',&
@@ -2247,6 +2342,7 @@
          epsij=eps(i,j)
          if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
            rrij=sigma(i,j)
          epsij=eps(i,j)
          if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
            rrij=sigma(i,j)
+!            print *,"SIGMA ZLA?",sigma(i,j)
           else
            rrij=rr0(i)+rr0(j)
           endif
           else
            rrij=rr0(i)+rr0(j)
           endif
@@ -2257,6 +2353,7 @@
          sigeps=dsign(1.0D0,epsij)
          epsij=dabs(epsij)
          aa_aq(i,j)=epsij*rrij*rrij
          sigeps=dsign(1.0D0,epsij)
          epsij=dabs(epsij)
          aa_aq(i,j)=epsij*rrij*rrij
+          print *,"ADASKO",epsij,rrij,expon
          bb_aq(i,j)=-sigeps*epsij*rrij
          aa_aq(j,i)=aa_aq(i,j)
          bb_aq(j,i)=bb_aq(i,j)
          bb_aq(i,j)=-sigeps*epsij*rrij
          aa_aq(j,i)=aa_aq(i,j)
          bb_aq(j,i)=bb_aq(i,j)
@@ -2268,7 +2365,7 @@
           aa_lip(j,i)=aa_lip(i,j)
           bb_lip(j,i)=bb_lip(i,j)
 !C          write(iout,*) aa_lip
           aa_lip(j,i)=aa_lip(i,j)
           bb_lip(j,i)=bb_lip(i,j)
 !C          write(iout,*) aa_lip
-         if (ipot.gt.2) then
+         if ((ipot.gt.2).and. (scelemode.eq.0)) then
            sigt1sq=sigma0(i)**2
            sigt2sq=sigma0(j)**2
            sigii1=sigii(i)
            sigt1sq=sigma0(i)**2
            sigt2sq=sigma0(j)**2
            sigii1=sigii(i)
@@ -2435,7 +2532,7 @@
         read(isidep_scbase,*) &
        dhead_scbasei(i,j), &
        dhead_scbasej(i,j), &
         read(isidep_scbase,*) &
        dhead_scbasei(i,j), &
        dhead_scbasej(i,j), &
-       rborn_scbasei(i,j),rborn_scbasej(j,i)
+       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,*) &
        (wdipdip_scbase(k,i,j),k=1,3), &
        (wqdip_scbase(k,i,j),k=1,2)
@@ -2500,6 +2597,77 @@
           aa_pepbase(j)=epsij*rrij*rrij
           bb_pepbase(j)=-sigeps*epsij*rrij
         enddo
           aa_pepbase(j)=epsij*rrij*rrij
           bb_pepbase(j)=-sigeps*epsij*rrij
         enddo
+!--------------READING SC PHOSPHATE------------------------------------- 
+      allocate(eps_scpho(ntyp_molec(1)))
+      allocate(sigma_scpho(ntyp_molec(1)))
+      allocate(chi_scpho(ntyp_molec(1),2))
+      allocate(chipp_scpho(ntyp_molec(1),2))
+      allocate(alphasur_scpho(4,ntyp_molec(1)))
+      allocate(sigmap1_scpho(ntyp_molec(1)))
+      allocate(sigmap2_scpho(ntyp_molec(1)))
+      allocate(chis_scpho(ntyp_molec(1),2))
+      allocate(wqq_scpho(ntyp_molec(1)))
+      allocate(wqdip_scpho(2,ntyp_molec(1)))
+      allocate(alphapol_scpho(ntyp_molec(1)))
+      allocate(epsintab_scpho(ntyp_molec(1)))
+      allocate(dhead_scphoi(ntyp_molec(1)))
+      allocate(rborn_scphoi(ntyp_molec(1)))
+      allocate(rborn_scphoj(ntyp_molec(1)))
+      allocate(alphi_scpho(ntyp_molec(1)))
+
+
+!      j=1
+       do j=1,ntyp_molec(1) ! without U then we will take T for U
+        write (*,*) "Im in scpho ", i, " ", j
+        read(isidep_scpho,*) &
+        eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
+        chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
+         write(*,*) "eps",eps_scpho(j)
+        read(isidep_scpho,*) &
+       (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
+       chis_scpho(j,1),chis_scpho(j,2)
+        read(isidep_scpho,*) &
+       (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
+        read(isidep_scpho,*) &
+         epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
+         alphi_scpho(j)
+       
+       END DO
+      allocate(aa_scpho(ntyp_molec(1)))
+      allocate(bb_scpho(ntyp_molec(1)))
+
+       do j=1,ntyp_molec(1)
+          epsij=eps_scpho(j)
+          rrij=sigma_scpho(j)
+!          r0(i,j)=rrij
+!          r0(j,i)=rrij
+          rrij=rrij**expon
+!          epsij=eps(i,j)
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_scpho(j)=epsij*rrij*rrij
+          bb_scpho(j)=-sigeps*epsij*rrij
+        enddo
+
+
+        read(isidep_peppho,*) &
+        eps_peppho,sigma_peppho
+        read(isidep_peppho,*) &
+       (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
+        read(isidep_peppho,*) &
+       (wqdip_peppho(k),k=1,2)
+
+          epsij=eps_peppho
+          rrij=sigma_peppho
+!          r0(i,j)=rrij
+!          r0(j,i)=rrij
+          rrij=rrij**expon
+!          epsij=eps(i,j)
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_peppho=epsij*rrij*rrij
+          bb_peppho=-sigeps*epsij*rrij
+
 
       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
       bad(:,:)=0.0D0
 
       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
       bad(:,:)=0.0D0
@@ -2603,6 +2771,31 @@
       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
         ' v3ss:',v3ss
       endif
       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
         ' v3ss:',v3ss
       endif
+      if (shield_mode.gt.0) then
+      pi=4.0D0*datan(1.0D0)
+!C VSolvSphere the volume of solving sphere
+      print *,pi,"pi"
+!C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
+!C there will be no distinction between proline peptide group and normal peptide
+!C group in case of shielding parameters
+      VSolvSphere=4.0/3.0*pi*(4.50d0)**3
+      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
+      write (iout,*) VSolvSphere,VSolvSphere_div
+!C long axis of side chain 
+      do i=1,ntyp
+      long_r_sidechain(i)=vbldsc0(1,i)
+!      if (scelemode.eq.0) then
+      short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
+      if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
+!      else
+!      short_r_sidechain(i)=sigma(i,i)
+!      endif
+      write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
+         sigma0(i) 
+      enddo
+      buff_shield=1.0d0
+      endif
+
       return
   111 write (iout,*) "Error reading bending energy parameters."
       goto 999
       return
   111 write (iout,*) "Error reading bending energy parameters."
       goto 999
@@ -2665,7 +2858,7 @@
       use control_data
       use compare_data
       use MPI_data
       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'
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.LOCAL'
@@ -2708,18 +2901,30 @@
 !      write (2,*) "UNRES_PDB",unres_pdb
       ires=0
       ires_old=0
 !      write (2,*) "UNRES_PDB",unres_pdb
       ires=0
       ires_old=0
+#ifdef WHAM_RUN
+      do i=1,nres
+       do j=1,5
+        itype(i,j)=0
+       enddo
+      enddo
+#endif
       nres=0
       iii=0
       lsecondary=.false.
       nhfrag=0
       nbfrag=0
       nres=0
       iii=0
       lsecondary=.false.
       nhfrag=0
       nbfrag=0
+      do j=1,5
+       nres_molec(j)=0
+      enddo
+      
+       
 !-----------------------------
       allocate(hfrag(2,maxres/3)) !(2,maxres/3)
       allocate(bfrag(4,maxres/3)) !(4,maxres/3)
 !-----------------------------
       allocate(hfrag(2,maxres/3)) !(2,maxres/3)
       allocate(bfrag(4,maxres/3)) !(4,maxres/3)
-
+      if(.not. allocated(istype)) allocate(istype(maxres))
       do i=1,100000
         read (ipdbin,'(a80)',end=10) card
       do i=1,100000
         read (ipdbin,'(a80)',end=10) card
-!       write (iout,'(a)') card
+       write (iout,'(a)') card
         if (card(:5).eq.'HELIX') then
           nhfrag=nhfrag+1
           lsecondary=.true.
         if (card(:5).eq.'HELIX') then
           nhfrag=nhfrag+1
           lsecondary=.true.
@@ -2829,7 +3034,7 @@
 !              nres_molec(molecule)=nres_molec(molecule)+1
             else
              molecule=2
 !              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)
 !              nres_molec(molecule)=nres_molec(molecule)+1
              read (card(19:19),'(a1)') sugar
              isugar=sugarcode(sugar,ires)
@@ -2899,11 +3104,19 @@
 !              print *,ires,istype(ires)
             endif
             if (unres_pdb) then
 !              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)
               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
               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. &
 !            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. &
@@ -2944,7 +3157,8 @@
            molecule=5
            nres_molec(molecule)=nres_molec(molecule)+1
            print *,"HERE",nres_molec(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)
+           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
            read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
           endif
         endif !atom
@@ -2954,9 +3168,10 @@
 ! Calculate dummy residue coordinates inside the "chain" of a multichain
 ! system
       nres=ires
 ! Calculate dummy residue coordinates inside the "chain" of a multichain
 ! system
       nres=ires
-      if ((ires_old.ne.ires).and.(molecule.ne.5)) &
+      if (((ires_old.ne.ires).and.(molecule.ne.5)) &
+        ) &
          nres_molec(molecule)=nres_molec(molecule)-2
          nres_molec(molecule)=nres_molec(molecule)-2
-!      print *,'I have', nres_molec(:)
+      print *,'I have',nres, nres_molec(:)
       
       do k=1,4 ! ions are without dummy 
        if (nres_molec(k).eq.0) cycle
       
       do k=1,4 ! ions are without dummy 
        if (nres_molec(k).eq.0) cycle
@@ -2988,47 +3203,47 @@
 ! 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
 ! 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'
 ! 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?'
 !            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
            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
          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
 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
-            call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
-            if (fail) then
-              e2(1)=0.0d0
-              e2(2)=1.0d0
-              e2(3)=0.0d0
-            endif
+!            call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
+!            if (fail) then
+!              e2(1)=0.0d0
+!              e2(2)=1.0d0
+!              e2(3)=0.0d0
+!            endif
             do j=1,3
               c(j,i)=c(j,i+1)-1.9d0*e2(j)
             enddo
             do j=1,3
               c(j,i)=c(j,i+1)-1.9d0*e2(j)
             enddo
-          else !unres_pdb
+!          else !unres_pdb
            do j=1,3
             dcj=(c(j,i+3)-c(j,i+2))/2.0
             if (dcj.eq.0) dcj=1.23591524223
             c(j,i)=c(j,i+1)-dcj
             c(j,nres+i)=c(j,i)
            enddo
            do j=1,3
             dcj=(c(j,i+3)-c(j,i+2))/2.0
             if (dcj.eq.0) dcj=1.23591524223
             c(j,i)=c(j,i+1)-dcj
             c(j,nres+i)=c(j,i)
            enddo
-          endif !unres_pdb
+!          endif !unres_pdb
          endif !itype(i+1,1).eq.ntyp1
         endif  !itype.eq.ntyp1
 
          endif !itype(i+1,1).eq.ntyp1
         endif  !itype.eq.ntyp1
 
@@ -3053,24 +3268,24 @@
           if (molecule.eq.5) molecule=molecprev
         itype(nres,molecule)=ntyp1_molec(molecule)
         nres_molec(molecule)=nres_molec(molecule)+1
           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
 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
-          call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
-          if (fail) then
-            e2(1)=0.0d0
-            e2(2)=1.0d0
-            e2(3)=0.0d0
-          endif
-          do j=1,3
-            c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
-          enddo
-        else
+!          call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
+!          if (fail) then
+!            e2(1)=0.0d0
+!            e2(2)=1.0d0
+!            e2(3)=0.0d0
+!          endif
+!          do j=1,3
+!            c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
+!          enddo
+!        else
         do j=1,3
           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
           c(j,nres)=c(j,nres-1)+dcj
           c(j,2*nres)=c(j,nres)
         enddo
         do j=1,3
           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
           c(j,nres)=c(j,nres-1)+dcj
           c(j,2*nres)=c(j,nres)
         enddo
-        endif
+!        endif
       endif
 !     print *,'I have',nres, nres_molec(:)
 
       endif
 !     print *,'I have',nres, nres_molec(:)
 
@@ -3172,10 +3387,10 @@
       if (lprn) then
       write (iout,'(/a)') &
         "Cartesian coordinates of the reference structure"
       if (lprn) then
       write (iout,'(/a)') &
         "Cartesian coordinates of the reference structure"
-      write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
+      write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
       do ires=1,nres
        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
       do ires=1,nres
-        write (iout,'(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
           (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
           (c(j,ires+nres),j=1,3)
       enddo
@@ -3186,7 +3401,7 @@
        write (iout,'(a)') &
          "Backbone and SC coordinates as read from the PDB"
        do ires=1,nres
        write (iout,'(a)') &
          "Backbone and SC coordinates as read from the PDB"
        do ires=1,nres
-        write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
+        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
           ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
           (c(j,nres+ires),j=1,3)
        enddo
@@ -3194,7 +3409,9 @@
 ! NOW LETS ROCK! SORTING
       allocate(c_temporary(3,2*nres))
       allocate(itype_temporary(nres,5))
 ! 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
       allocate(istype_temp(nres))
        itype_temporary(:,:)=0
       seqalingbegin=1
@@ -3225,36 +3442,36 @@
          istype(i)=istype_temp(i)
         enddo
        enddo
          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
 ! 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"
 
       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
        "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
           (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
           (c(j,ires+nres),j=1,3)
       enddo
@@ -3342,7 +3559,7 @@
         write (iout,*) "symetr", symetr
       do i=1,nres
       lll=lll+1
         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
       if (i.gt.1) then
       if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
       chain_length=lll-1
@@ -3371,7 +3588,7 @@
 !       write (iout,*) "spraw lancuchy",chain_length,symetr
 !       do i=1,4
 !         do kkk=1,chain_length
 !       write (iout,*) "spraw lancuchy",chain_length,symetr
 !       do i=1,4
 !         do kkk=1,chain_length
-!           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
+!           write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
 !         enddo
 !        enddo
 ! enddiagnostic
 !         enddo
 !        enddo
 ! enddiagnostic
@@ -3394,7 +3611,7 @@
         cou=0
         do kkk=1,symetr
          icha=tabperm(i,kkk)
         cou=0
         do kkk=1,symetr
          icha=tabperm(i,kkk)
-!         write (iout,*) i,icha
+         write (iout,*) i,icha
          do lll=1,chain_length
           cou=cou+1
            if (cou.le.nres) then
          do lll=1,chain_length
           cou=cou+1
            if (cou.le.nres) then
@@ -3402,7 +3619,7 @@
             kupa=mod(lll,chain_length)
             iprzes=(kkk-1)*chain_length+lll
             if (kupa.eq.0) kupa=chain_length
             kupa=mod(lll,chain_length)
             iprzes=(kkk-1)*chain_length+lll
             if (kupa.eq.0) kupa=chain_length
-!            write (iout,*) "kupa", kupa
+            write (iout,*) "kupa", kupa
             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
           enddo
             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
           enddo
@@ -3421,11 +3638,11 @@
       cref(3,i,kkk),cref(1,nres+i,kkk),&
       cref(2,nres+i,kkk),cref(3,nres+i,kkk)
       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')
                 '     centroid coordinates'/ &
                 '       ', 6X,'X',11X,'Y',11X,'Z', &
                                 10X,'X',11X,'Y',11X,'Z')
-  110 format (a,'(',i3,')',6f12.5)
+  110 format (a,'(',i5,')',6f12.5)
      
       enddo
 !c enddiag
      
       enddo
 !c enddiag
@@ -3555,6 +3772,13 @@
       call reada(controlcard,'BOXY',boxysize,100.0d0)
       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
       call readi(controlcard,'TUBEMOD',tubemode,0)
       call reada(controlcard,'BOXY',boxysize,100.0d0)
       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
       call readi(controlcard,'TUBEMOD',tubemode,0)
+      print *,"SCELE",scelemode
+      call readi(controlcard,"SCELEMODE",scelemode,0)
+      print *,"SCELE",scelemode
+
+! elemode = 0 is orignal UNRES electrostatics
+! elemode = 1 is "Momo" potentials in progress
+! elemode = 2 is in development EVALD
       write (iout,*) TUBEmode,"TUBEMODE"
       if (TUBEmode.gt.0) then
        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
       write (iout,*) TUBEmode,"TUBEMODE"
       if (TUBEmode.gt.0) then
        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
@@ -3692,22 +3916,22 @@
        write (iout,'(2a)') diagmeth(kdiag),&
         ' routine used to diagonalize matrices.'
       if (shield_mode.gt.0) then
        write (iout,'(2a)') diagmeth(kdiag),&
         ' routine used to diagonalize matrices.'
       if (shield_mode.gt.0) then
-      pi=3.141592d0
+       pi=4.0D0*datan(1.0D0)
 !C VSolvSphere the volume of solving sphere
 !C VSolvSphere the volume of solving sphere
-!C      print *,pi,"pi"
+      print *,pi,"pi"
 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
 !C there will be no distinction between proline peptide group and normal peptide
 !C group in case of shielding parameters
 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
 !C there will be no distinction between proline peptide group and normal peptide
 !C group in case of shielding parameters
-      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
-      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
+      VSolvSphere=4.0/3.0*pi*(4.50d0)**3
+      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
       write (iout,*) VSolvSphere,VSolvSphere_div
 !C long axis of side chain 
       write (iout,*) VSolvSphere,VSolvSphere_div
 !C long axis of side chain 
-      do i=1,ntyp
-      long_r_sidechain(i)=vbldsc0(1,i)
-      short_r_sidechain(i)=sigma0(i)
-      write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
-         sigma0(i) 
-      enddo
+!      do i=1,ntyp
+!      long_r_sidechain(i)=vbldsc0(1,i)
+!      short_r_sidechain(i)=sigma0(i)
+!      write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
+!         sigma0(i) 
+!      enddo
       buff_shield=1.0d0
       endif
       return
       buff_shield=1.0d0
       endif
       return
@@ -3867,6 +4091,10 @@
       large = index(controlcard,"LARGE").gt.0
       print_compon = index(controlcard,"PRINT_COMPON").gt.0
       rattle = index(controlcard,"RATTLE").gt.0
       large = index(controlcard,"LARGE").gt.0
       print_compon = index(controlcard,"PRINT_COMPON").gt.0
       rattle = index(controlcard,"RATTLE").gt.0
+      preminim=(index(controlcard,'PREMINIM').gt.0)
+      write (iout,*) "PREMINIM ",preminim
+      dccart=(index(controlcard,'CART').gt.0)
+      if (preminim) call read_minim
 !  if performing umbrella sampling, fragments constrained are read from the fragment file 
       nset=0
       if(usampl) then
 !  if performing umbrella sampling, fragments constrained are read from the fragment file 
       nset=0
       if(usampl) then
@@ -4503,6 +4731,11 @@
       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')
       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('LIPTRANPAR',liptranname)
       open (iliptranpar,file=liptranname,status='old',action='read')
UNRESPACK 4.x