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 e6b7f0a..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,11 +736,11 @@
       character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
       logical :: lprint,LaTeX
       real(kind=8),dimension(3,3,maxlob) :: blower     !(3,3,maxlob)
       character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
       logical :: lprint,LaTeX
       real(kind=8),dimension(3,3,maxlob) :: blower     !(3,3,maxlob)
-      real(kind=8),dimension(13) :: b
+      real(kind=8),dimension(13) :: buse
       character(len=3) :: lancuch      !,ucase
 !el  local variables
       integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
       character(len=3) :: lancuch      !,ucase
 !el  local variables
       integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
-      integer :: maxinter,junk,kk,ii
+      integer :: maxinter,junk,kk,ii,ncatprotparm
       real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
                 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
                 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
       real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
                 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
                 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
@@ -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),&
@@ -854,6 +865,19 @@
           enddo
         enddo
       endif
           enddo
         enddo
       endif
+            do i=1,ntyp_molec(5)
+             read(iion,*) msc(i,5),restok(i,5)
+             print *,msc(i,5),restok(i,5)
+            enddo
+            ip(5)=0.2
+!            isc(5)=0.2
+            read (iion,*) ncatprotparm
+            allocate(catprm(ncatprotparm,4))
+            do k=1,4
+            read (iion,*)  (catprm(i,k),i=1,ncatprotparm)
+            enddo
+            print *, catprm
+!            read (iion,*) (vcatprm(k),k=1,ncatprotpram)
 !----------------------------------------------------
       allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
       allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp))     !(-ntyp:ntyp)
 !----------------------------------------------------
       allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
       allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp))     !(-ntyp:ntyp)
@@ -1412,6 +1436,24 @@
          enddo  
        endif
       enddo
          enddo  
        endif
       enddo
+!---------reading nucleic acid parameters for rotamers-------------------
+      allocate(sc_parmin_nucl(9,ntyp_molec(2)))      !(maxsccoef,ntyp)
+      do i=1,ntyp_molec(2)
+        read (irotam_nucl,*,end=112,err=112)
+        do j=1,9
+          read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
+        enddo
+      enddo
+      close(irotam_nucl)
+      if (lprint) then
+        write (iout,*)
+        write (iout,*) "Base rotamer parameters"
+        do i=1,ntyp_molec(2)
+          write (iout,'(a)') restyp(i,2)
+          write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
+        enddo
+      endif
+
 !
 ! Read the parameters of the probability distribution/energy expression
 ! of the side chains.
 !
 ! Read the parameters of the probability distribution/energy expression
 ! of the side chains.
@@ -1866,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
@@ -1910,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
@@ -2115,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
@@ -2142,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)
@@ -2170,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),&
@@ -2181,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
@@ -2199,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)
@@ -2207,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   ',&
@@ -2216,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
@@ -2226,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)
@@ -2237,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)
@@ -2372,6 +2500,174 @@
        sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
       write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
       enddo
        sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
       write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
       enddo
+!-----------------READING SC BASE POTENTIALS-----------------------------
+      allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))      
+      allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
+      allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
+      allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
+      allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
+      allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
+      allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
+      allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
+      allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
+      allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
+      allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
+      allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
+      allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
+      allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
+      allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
+      allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
+
+
+      do i=1,ntyp_molec(1)
+       do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
+        write (*,*) "Im in ", i, " ", j
+        read(isidep_scbase,*) &
+        eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
+        chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
+         write(*,*) "eps",eps_scbase(i,j)
+        read(isidep_scbase,*) &
+       (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
+       chis_scbase(i,j,1),chis_scbase(i,j,2)
+        read(isidep_scbase,*) &
+       dhead_scbasei(i,j), &
+       dhead_scbasej(i,j), &
+       rborn_scbasei(i,j),rborn_scbasej(i,j)
+        read(isidep_scbase,*) &
+       (wdipdip_scbase(k,i,j),k=1,3), &
+       (wqdip_scbase(k,i,j),k=1,2)
+        read(isidep_scbase,*) &
+       alphapol_scbase(i,j), &
+       epsintab_scbase(i,j) 
+       END DO
+      END DO
+      allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
+      allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
+
+      do i=1,ntyp_molec(1)
+       do j=1,ntyp_molec(2)-1 
+          epsij=eps_scbase(i,j)
+          rrij=sigma_scbase(i,j)
+!          r0(i,j)=rrij
+!          r0(j,i)=rrij
+          rrij=rrij**expon
+!          epsij=eps(i,j)
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_scbase(i,j)=epsij*rrij*rrij
+          bb_scbase(i,j)=-sigeps*epsij*rrij
+        enddo
+       enddo
+!-----------------READING PEP BASE POTENTIALS-------------------
+      allocate(eps_pepbase(ntyp_molec(2)))
+      allocate(sigma_pepbase(ntyp_molec(2)))
+      allocate(chi_pepbase(ntyp_molec(2),2))
+      allocate(chipp_pepbase(ntyp_molec(2),2))
+      allocate(alphasur_pepbase(4,ntyp_molec(2)))
+      allocate(sigmap1_pepbase(ntyp_molec(2)))
+      allocate(sigmap2_pepbase(ntyp_molec(2)))
+      allocate(chis_pepbase(ntyp_molec(2),2))
+      allocate(wdipdip_pepbase(3,ntyp_molec(2)))
+
+
+       do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
+        write (*,*) "Im in ", i, " ", j
+        read(isidep_pepbase,*) &
+        eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
+        chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
+         write(*,*) "eps",eps_pepbase(j)
+        read(isidep_pepbase,*) &
+       (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
+       chis_pepbase(j,1),chis_pepbase(j,2)
+        read(isidep_pepbase,*) &
+       (wdipdip_pepbase(k,j),k=1,3)
+       END DO
+      allocate(aa_pepbase(ntyp_molec(2)))
+      allocate(bb_pepbase(ntyp_molec(2)))
+
+       do j=1,ntyp_molec(2)-1
+          epsij=eps_pepbase(j)
+          rrij=sigma_pepbase(j)
+!          r0(i,j)=rrij
+!          r0(j,i)=rrij
+          rrij=rrij**expon
+!          epsij=eps(i,j)
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_pepbase(j)=epsij*rrij*rrij
+          bb_pepbase(j)=-sigeps*epsij*rrij
+        enddo
+!--------------READING SC PHOSPHATE------------------------------------- 
+      allocate(eps_scpho(ntyp_molec(1)))
+      allocate(sigma_scpho(ntyp_molec(1)))
+      allocate(chi_scpho(ntyp_molec(1),2))
+      allocate(chipp_scpho(ntyp_molec(1),2))
+      allocate(alphasur_scpho(4,ntyp_molec(1)))
+      allocate(sigmap1_scpho(ntyp_molec(1)))
+      allocate(sigmap2_scpho(ntyp_molec(1)))
+      allocate(chis_scpho(ntyp_molec(1),2))
+      allocate(wqq_scpho(ntyp_molec(1)))
+      allocate(wqdip_scpho(2,ntyp_molec(1)))
+      allocate(alphapol_scpho(ntyp_molec(1)))
+      allocate(epsintab_scpho(ntyp_molec(1)))
+      allocate(dhead_scphoi(ntyp_molec(1)))
+      allocate(rborn_scphoi(ntyp_molec(1)))
+      allocate(rborn_scphoj(ntyp_molec(1)))
+      allocate(alphi_scpho(ntyp_molec(1)))
+
+
+!      j=1
+       do j=1,ntyp_molec(1) ! without U then we will take T for U
+        write (*,*) "Im in scpho ", i, " ", j
+        read(isidep_scpho,*) &
+        eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
+        chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
+         write(*,*) "eps",eps_scpho(j)
+        read(isidep_scpho,*) &
+       (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
+       chis_scpho(j,1),chis_scpho(j,2)
+        read(isidep_scpho,*) &
+       (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
+        read(isidep_scpho,*) &
+         epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
+         alphi_scpho(j)
+       
+       END DO
+      allocate(aa_scpho(ntyp_molec(1)))
+      allocate(bb_scpho(ntyp_molec(1)))
+
+       do j=1,ntyp_molec(1)
+          epsij=eps_scpho(j)
+          rrij=sigma_scpho(j)
+!          r0(i,j)=rrij
+!          r0(j,i)=rrij
+          rrij=rrij**expon
+!          epsij=eps(i,j)
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_scpho(j)=epsij*rrij*rrij
+          bb_scpho(j)=-sigeps*epsij*rrij
+        enddo
+
+
+        read(isidep_peppho,*) &
+        eps_peppho,sigma_peppho
+        read(isidep_peppho,*) &
+       (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
+        read(isidep_peppho,*) &
+       (wqdip_peppho(k),k=1,2)
+
+          epsij=eps_peppho
+          rrij=sigma_peppho
+!          r0(i,j)=rrij
+!          r0(j,i)=rrij
+          rrij=rrij**expon
+!          epsij=eps(i,j)
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_peppho=epsij*rrij*rrij
+          bb_peppho=-sigeps*epsij*rrij
+
 
       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
       bad(:,:)=0.0D0
 
       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
       bad(:,:)=0.0D0
@@ -2475,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
@@ -2537,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'
@@ -2560,7 +2881,7 @@
       character(len=80) :: card
       real(kind=8),dimension(3,20) :: sccor
       integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin      !rescode,
       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
       character*1 :: sugar
       real(kind=8) :: cou
       real(kind=8),dimension(3) :: ccc
@@ -2580,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.
@@ -2697,10 +3030,11 @@
              if (atom.eq.'CA  '.or.atom.eq.'N   ') then
              molecule=1
               itype(ires,molecule)=rescode(ires,res,0,molecule)
              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
 !              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)
@@ -2749,9 +3083,11 @@
             do j=1,3
               c(j,ires)=c(j,ires)+ccc(j)/5.0
             enddo
             do j=1,3
               c(j,ires)=c(j,ires)+ccc(j)/5.0
             enddo
+             print *,counter,molecule
              if (counter.eq.5) then
 !            iii=iii+1
               nres_molec(molecule)=nres_molec(molecule)+1
              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
 !            do j=1,3
 !              sccor(j,iii)=c(j,ires)
 !            enddo
@@ -2768,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. &
@@ -2791,7 +3135,16 @@
               endif
           endif
         else if ((ions).and.(card(1:6).eq.'HETATM')) then
               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
           read (card(12:16),*) atom
           print *,"HETATOM", atom
           read (card(18:20),'(a3)') res
@@ -2803,7 +3156,9 @@
            if (molecule.ne.5) molecprev=molecule
            molecule=5
            nres_molec(molecule)=nres_molec(molecule)+1
            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
            read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
           endif
         endif !atom
@@ -2813,8 +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) 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
       
       do k=1,4 ! ions are without dummy 
        if (nres_molec(k).eq.0) cycle
@@ -2846,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
 
@@ -2911,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(:)
 
@@ -3030,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
@@ -3044,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
@@ -3052,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
@@ -3065,6 +3424,7 @@
 
           enddo
           itype_temporary(seqalingbegin,k)=itype(i,k)
 
           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_temp(seqalingbegin)=istype(i)
           molnum(seqalingbegin)=k
           seqalingbegin=seqalingbegin+1
@@ -3082,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
@@ -3199,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
@@ -3228,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
@@ -3251,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
@@ -3259,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
@@ -3278,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
@@ -3412,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)
@@ -3549,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
@@ -3689,7 +4056,7 @@
 !      character(len=80) :: ucase
       character(len=320) :: controlcard
 !el local variables
 !      character(len=80) :: ucase
       character(len=320) :: controlcard
 !el local variables
-      integer :: i
+      integer :: i,j
       real(kind=8) :: eta
 
       call card_concat(controlcard,.true.)
       real(kind=8) :: eta
 
       call card_concat(controlcard,.true.)
@@ -3724,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
@@ -3802,21 +4173,27 @@
         if(me.eq.king.or..not.out1file) &
          write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
           eta
         if(me.eq.king.or..not.out1file) &
          write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
           eta
-        gamp=scal_fric*(pstok(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
         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 
+        enddo
+
         if(me.eq.king.or..not.out1file)then
          write (iout,'(/2a/)') &
          "Radii of site types and friction coefficients and std's of",&
          " stochastic forces of fully exposed sites"
         if(me.eq.king.or..not.out1file)then
          write (iout,'(/2a/)') &
          "Radii of site types and friction coefficients and std's of",&
          " stochastic forces of fully exposed sites"
-         write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
+         write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
          do i=1,ntyp
           write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,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
          enddo
         endif
       else if (tbf) then
@@ -4259,6 +4636,8 @@
       open (iliptranpar,file=liptranname,status='old',action='read')
       call getenv_loc('TUBEPAR',tubename)
       open (itube,file=tubename,status='old',action='read')
       open (iliptranpar,file=liptranname,status='old',action='read')
       call getenv_loc('TUBEPAR',tubename)
       open (itube,file=tubename,status='old',action='read')
+      call getenv_loc('IONPAR',ionname)
+      open (iion,file=ionname,status='old',action='read')
 
 !      print *,"Processor",myrank," opened file ISIDEP" 
 !      print *,"Processor",myrank," opened parameter files" 
 
 !      print *,"Processor",myrank," opened file ISIDEP" 
 !      print *,"Processor",myrank," opened parameter files" 
@@ -4303,6 +4682,8 @@
       open (iliptranpar,file=liptranname,status='old')
       call getenv_loc('TUBEPAR',tubename)
       open (itube,file=tubename,status='old')
       open (iliptranpar,file=liptranname,status='old')
       call getenv_loc('TUBEPAR',tubename)
       open (itube,file=tubename,status='old')
+      call getenv_loc('IONPAR',ionname)
+      open (iion,file=ionname,status='old')
 #else
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
         readonly)
 #else
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
         readonly)
@@ -4346,11 +4727,22 @@
       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')
       open (itordp_nucl,file=tordname_nucl,status='old',action='read')
       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
       open (isidep_nucl,file=sidename_nucl,status='old',action='read')
+      call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
+      open (isidep_scbase,file=sidename_scbase,status='old',action='read')
+      call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
+      open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
+      call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
+      open (isidep_scpho,file=pepname_scpho,status='old',action='read')
+      call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
+      open (isidep_peppho,file=pepname_peppho,status='old',action='read')
+
 
       call getenv_loc('LIPTRANPAR',liptranname)
       open (iliptranpar,file=liptranname,status='old',action='read')
       call getenv_loc('TUBEPAR',tubename)
       open (itube,file=tubename,status='old',action='read')
 
       call getenv_loc('LIPTRANPAR',liptranname)
       open (iliptranpar,file=liptranname,status='old',action='read')
       call getenv_loc('TUBEPAR',tubename)
       open (itube,file=tubename,status='old',action='read')
+      call getenv_loc('IONPAR',ionname)
+      open (iion,file=ionname,status='old',action='read')
 
 #ifndef CRYST_SC
       call getenv_loc('ROTPARPDB',rotname_pdb)
 
 #ifndef CRYST_SC
       call getenv_loc('ROTPARPDB',rotname_pdb)
UNRESPACK 4.x