reading from extended mutliple sequence
[unres4.git] / source / unres / io_config.f90
index da0414d..277b6ba 100644 (file)
 
 !     write (iout,100)
 !      do i=1,nres
-!        write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
+!        write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
 !          c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
 !      enddo
 !  100 format (//'              alpha-carbon coordinates       ',&
       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,&
-                res1
+                res1,epsijlip,epspeptube,epssctube,sigmapeptube,      &
+                sigmasctube
       integer :: ichir1,ichir2
 !      real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
 !el      allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
       allocate(isc(ntyp+1)) !(ntyp+1)
       allocate(restok(ntyp+1)) !(ntyp+1)
       allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
-
+      allocate(long_r_sidechain(ntyp))
+      allocate(short_r_sidechain(ntyp))
       dsc(:)=0.0d0
       dsc_inv(:)=0.0d0
 
         endif
       enddo
 #else
-      read (ibond,*) junk,vbldpDUM,vbldp0,akp,rjunk,mp,ip,pstok
+      read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp,ip,pstok
       do i=1,ntyp
         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
          j=1,nbondterm(i)),msc(i),isc(i),restok(i)
          'inertia','Pstok'
         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
         do i=1,ntyp
-          write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
+          write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
             vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
           do j=2,nbondterm(i)
             write (iout,'(13x,3f10.5)') &
       theta0(:)=0.0D0
       sig0(:)=0.0D0
       sigc0(:)=0.0D0
+      allocate(liptranene(ntyp))
+!C reading lipid parameters
+      write (iout,*) "iliptranpar",iliptranpar
+      call flush(iout)
+       read(iliptranpar,*) pepliptran
+       print *,pepliptran
+       do i=1,ntyp
+       read(iliptranpar,*) liptranene(i)
+        print *,liptranene(i)
+       enddo
+       close(iliptranpar)
 
 #ifdef CRYST_THETA
 !
        '    ATHETA0   ','         A1   ','        A2    ',&
        '        B1    ','         B2   '        
         do i=1,ntyp
-          write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
+          write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
               a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
         enddo
         write (iout,'(/a/9x,5a/79(1h-))') &
        '     ALPH0    ','      ALPH1   ','     ALPH2    ',&
        '     ALPH3    ','    SIGMA0C   '        
         do i=1,ntyp
-          write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
+          write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
             (polthet(j,i),j=0,3),sigc0(i) 
         enddo
         write (iout,'(/a/9x,5a/79(1h-))') &
        '    THETA0    ','     SIGMA0   ','        G1    ',&
        '        G2    ','         G3   '        
         do i=1,ntyp
-          write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),&
+          write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
              sig0(i),(gthet(j,i),j=1,3)
         enddo
        else
        '     theta0   ','    a1*10^2   ','   a2*10^2    ',&
        '   b1*10^1    ','    b2*10^1   '        
         do i=1,ntyp
-          write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
+          write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
               a0thet(i),(100*athet(j,i,1,1),j=1,2),&
               (10*bthet(j,i,1,1),j=1,2)
         enddo
        ' alpha0       ','  alph1       ',' alph2        ',&
        ' alhp3        ','   sigma0c    '        
        do i=1,ntyp
-          write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
+          write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
             (polthet(j,i),j=0,3),sigc0(i) 
        enddo
        write (iout,'(/a/9x,5a/79(1h-))') &
        '    theta0    ','  sigma0*10^2 ','      G1*10^-1',&
        '        G2    ','   G3*10^1    '        
        do i=1,ntyp
-          write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
+          write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
              100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
        enddo
       endif
          nlobi=nlob(i)
           if (nlobi.gt.0) then
             if (LaTeX) then
-              write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
+              write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
                ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
                                    'log h',(bsc(j,i),j=1,nlobi)
       do i=-ntyp,-1
        itortyp(i)=-itortyp(-i)
       enddo
-!      itortyp(ntyp1)=ntortyp
-!      itortyp(-ntyp1)=-ntortyp
+      itortyp(ntyp1)=ntortyp
+      itortyp(-ntyp1)=-ntortyp
       do iblock=1,2 
       write (iout,*) 'ntortyp',ntortyp
       do i=0,ntortyp-1
       allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
       allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
       allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
-
+      allocate(epslip(ntyp,ntyp))
       augm(:,:)=0.0D0
       chip(:)=0.0D0
       alp(:)=0.0D0
          call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
          write (iout,'(/a)') 'One-body parameters:'
          write (iout,'(a,4x,a)') 'residue','sigma'
-         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
+         write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
         endif
 !      goto 50
 !----------------------- LJK potential --------------------------------
          call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
          write (iout,'(/a)') 'One-body parameters:'
          write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
-          write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
+          write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
                 i=1,ntyp)
         endif
 !      goto 50
         read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
         read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
         read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
+        do i=1,ntyp
+         read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
+        enddo
+
 ! For the GB potential convert sigma'**2 into chi'
         if (ipot.eq.4) then
          do i=1,ntyp
          write (iout,'(/a)') 'One-body parameters:'
          write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',&
                '    chip  ','    alph  '
-         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
+         write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
                              chip(i),alp(i),i=1,ntyp)
         endif
 !      goto 50
           write (iout,'(/a)') 'One-body parameters:'
           write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',&
               's||/s_|_^2','    chip  ','    alph  '
-          write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
+          write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
                    sigii(i),chip(i),alp(i),i=1,ntyp)
         endif
        case default
 ! Calculate the "working" parameters of SC interactions.
 
 !el from module energy - COMMON.INTERACT-------
-      allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
+      allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
+      allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
       allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
-      aa(:,:)=0.0D0
-      bb(:,:)=0.0D0
+      allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
+        dcavtub(ntyp1))
+      allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
+        tubetranene(ntyp1))
+      aa_aq(:,:)=0.0D0
+      bb_aq(:,:)=0.0D0
+      aa_lip(:,:)=0.0D0
+      bb_lip(:,:)=0.0D0
       chi(:,:)=0.0D0
       sigma(:,:)=0.0D0
       r0(:,:)=0.0D0
+      acavtub(:)=0.0d0
+      bcavtub(:)=0.0d0
+      ccavtub(:)=0.0d0
+      dcavtub(:)=0.0d0
+      sc_aa_tube_par(:)=0.0d0
+      sc_bb_tube_par(:)=0.0d0
+
 !--------------------------------
 
       do i=2,ntyp
         do j=1,i-1
           eps(i,j)=eps(j,i)
+          epslip(i,j)=epslip(j,i)
         enddo
       enddo
       do i=1,ntyp
          epsij=eps(i,j)
          sigeps=dsign(1.0D0,epsij)
          epsij=dabs(epsij)
-         aa(i,j)=epsij*rrij*rrij
-         bb(i,j)=-sigeps*epsij*rrij
-         aa(j,i)=aa(i,j)
-         bb(j,i)=bb(i,j)
+         aa_aq(i,j)=epsij*rrij*rrij
+         bb_aq(i,j)=-sigeps*epsij*rrij
+         aa_aq(j,i)=aa_aq(i,j)
+         bb_aq(j,i)=bb_aq(i,j)
+          epsijlip=epslip(i,j)
+          sigeps=dsign(1.0D0,epsijlip)
+          epsijlip=dabs(epsijlip)
+          aa_lip(i,j)=epsijlip*rrij*rrij
+          bb_lip(i,j)=-sigeps*epsijlip*rrij
+          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
            sigt1sq=sigma0(i)**2
            sigt2sq=sigma0(j)**2
           endif
          if (lprint) then
             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
-            restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),&
+            restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
             sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
          endif
         enddo
       enddo
+      write(iout,*) "tube param"
+      read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
+      ccavtubpep,dcavtubpep,tubetranenepep
+      sigmapeptube=sigmapeptube**6
+      sigeps=dsign(1.0D0,epspeptube)
+      epspeptube=dabs(epspeptube)
+      pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
+      pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
+      write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
+      do i=1,ntyp
+       read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
+      ccavtub(i),dcavtub(i),tubetranene(i)
+       sigmasctube=sigmasctube**6
+       sigeps=dsign(1.0D0,epssctube)
+       epssctube=dabs(epssctube)
+       sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
+       sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
+      write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
+      enddo
 
       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
       bad(:,:)=0.0D0
 ! End current chain
           ires_old=ires+2
           ishift1=ishift1+1
-          itype(ires_old)=ntyp1
-          itype(ires_old-1)=ntyp1
+          itype(ires_old,1)=ntyp1
+          itype(ires_old-1,1)=ntyp1
           ibeg=2
 !          write (iout,*) "Chain ended",ires,ishift,ires_old
           if (unres_pdb) then
             enddo
           else
             call sccenter(ires,iii,sccor)
+!          iii=0
           endif
           iii=0
         endif
 !              write (iout,*) "Calculating sidechain center iii",iii
               if (unres_pdb) then
                 do j=1,3
-                  dc(j,ires+nres)=sccor(j,iii)
+                  dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
                 enddo
               else
                 call sccenter(ires_old,iii,sccor)
               ishift=ires-1
               if (res.ne.'GLY' .and. res.ne. 'ACE') then
                 ishift=ishift-1
-                itype(1)=ntyp1
+                itype(1,1)=ntyp1
               endif
               ires=ires-ishift+ishift1
               ires_old=ires
               ires_old=ires
             endif
             if (res.eq.'ACE' .or. res.eq.'NHE') then
-              itype(ires)=10
+              itype(ires,1)=10
             else
-              itype(ires)=rescode(ires,res,0)
+              itype(ires,1)=rescode(ires,res,0,1)
             endif
           else
             ires=ires-ishift+ishift1
 !            write (iout,*) "backbone ",atom
 #ifdef DEBUG
             write (iout,'(2i3,2x,a,3f8.3)') &
-            ires,itype(ires),res,(c(j,ires),j=1,3)
+            ires,itype(ires,1),res,(c(j,ires),j=1,3)
 #endif
             iii=iii+1
             do j=1,3
               sccor(j,iii)=c(j,ires)
             enddo
-!            write (*,*) card(23:27),ires,itype(ires)
+!            write (*,*) card(23:27),ires,itype(ires,1)
           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
                    atom.ne.'N' .and. atom.ne.'C' .and. &
                    atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
 ! system
       nres=ires
       do i=2,nres-1
-!        write (iout,*) i,itype(i)
-!        if (itype(i).eq.ntyp1) then
-!          write (iout,*) "dummy",i,itype(i)
+!        write (iout,*) i,itype(i,1)
+!        if (itype(i,1).eq.ntyp1) then
+!          write (iout,*) "dummy",i,itype(i,1)
 !          do j=1,3
 !            c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
 !            c(j,i)=(c(j,i-1)+c(j,i+1))/2
 !            dc(j,i)=c(j,i)
 !          enddo
 !        endif
-        if (itype(i).eq.ntyp1) then
-         if (itype(i+1).eq.ntyp1) then
+        if (itype(i,1).eq.ntyp1) then
+         if (itype(i+1,1).eq.ntyp1) then
 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
-! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
-! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
+! 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
 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
 !            print *,i,'tu dochodze'
              c(j,nres+i)=c(j,i)
            enddo
           endif   !unres_pdb
-         else     !itype(i+1).eq.ntyp1
+         else     !itype(i+1,1).eq.ntyp1
           if (unres_pdb) then
 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
             call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
             c(j,nres+i)=c(j,i)
            enddo
           endif !unres_pdb
-         endif !itype(i+1).eq.ntyp1
+         endif !itype(i+1,1).eq.ntyp1
         endif  !itype.eq.ntyp1
 
       enddo
 !      nres=ires
       nsup=nres
       nstart_sup=1
-      if (itype(nres).ne.10) then
+      if (itype(nres,1).ne.10) then
         nres=nres+1
-        itype(nres)=ntyp1
+        itype(nres,1)=ntyp1
         if (unres_pdb) then
 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
             e2(3)=0.0d0
           endif
           do j=1,3
-            c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
+            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)
+          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
         c(j,nres+1)=c(j,1)
         c(j,2*nres)=c(j,nres)
       enddo
-      if (itype(1).eq.ntyp1) then
+      if (itype(1,1).eq.ntyp1) then
         nsup=nsup-1
         nstart_sup=2
         if (unres_pdb) then
             e2(3)=0.0d0
           endif
           do j=1,3
-            c(j,1)=c(j,2)-3.8d0*e2(j)
+            c(j,1)=c(j,2)-1.9d0*e2(j)
           enddo
         else
         do j=1,3
-          dcj=c(j,4)-c(j,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
        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
       do ires=1,nres
         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
-          restyp(itype(ires)),ires,(c(j,ires),j=1,3),&
+          restyp(itype(ires,1),1),ires,(c(j,ires),j=1,3),&
           (c(j,ires+nres),j=1,3)
       enddo
       endif
          "Backbone and SC coordinates as read from the PDB"
        do ires=1,nres
         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
-          ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),&
+          ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
           (c(j,nres+ires),j=1,3)
        enddo
       endif
       lll=lll+1
 !c      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
       if (i.gt.1) then
-      if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
+      if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
       chain_length=lll-1
       kkk=kkk+1
 !       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
 ! enddiagnostic
 ! makes copy of chains
         write (iout,*) "symetr", symetr
+      do j=1,3
+      dc(j,0)=c(j,1)
+      enddo
 
       if (symetr.gt.1) then
        call permut(symetr)
       do kkk=1,nperm
       write (iout,*) "nowa struktura", nperm
       do i=1,nres
-        write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),&
+        write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
       cref(2,i,kkk),&
       cref(3,i,kkk),cref(1,nres+i,kkk),&
       cref(2,nres+i,kkk),cref(3,nres+i,kkk)
       character(len=640) :: controlcard
 
       real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
-                 
+      integer i                 
 
       nglob_csa=0
       eglob_csa=1d99
       timem=timlim
       modecalc=0
       call reada(controlcard,"T_BATH",t_bath,300.0d0)
+!C SHIELD keyword sets if the shielding effect of side-chains is used
+!C 0 denots no shielding is used all peptide are equally despite the 
+!C solvent accesible area
+!C 1 the newly introduced function
+!C 2 reseved for further possible developement
+      call readi(controlcard,'SHIELD',shield_mode,0)
+!C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+        write(iout,*) "shield_mode",shield_mode
 !C  Varibles set size of box
+      with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
+      protein=index(controlcard,"PROTEIN").gt.0
+      ions=index(controlcard,"IONS").gt.0
+      nucleic=index(controlcard,"NUCLEIC").gt.0
+      write (iout,*) "with_theta_constr ",with_theta_constr
+      AFMlog=(index(controlcard,'AFM'))
+      selfguide=(index(controlcard,'SELFGUIDE'))
+      print *,'AFMlog',AFMlog,selfguide,"KUPA"
+      call readi(controlcard,'GENCONSTR',genconstr,0)
       call reada(controlcard,'BOXX',boxxsize,100.0d0)
       call reada(controlcard,'BOXY',boxysize,100.0d0)
       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+      call readi(controlcard,'TUBEMOD',tubemode,0)
+      write (iout,*) TUBEmode,"TUBEMODE"
+      if (TUBEmode.gt.0) then
+       call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
+       call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
+       call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
+       call reada(controlcard,"RTUBE",tubeR0,0.0d0)
+       call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
+       call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
+       call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
+       buftubebot=bordtubebot+tubebufthick
+       buftubetop=bordtubetop-tubebufthick
+      endif
+
 ! CUTOFFF ON ELECTROSTATICS
       call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
+      write(iout,*) "R_CUT_ELE=",r_cut_ele
+! Lipidic parameters
+      call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
+      call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
+      if (lipthick.gt.0.0d0) then
+       bordliptop=(boxzsize+lipthick)/2.0
+       bordlipbot=bordliptop-lipthick
+      if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
+      write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
+      buflipbot=bordlipbot+lipbufthick
+      bufliptop=bordliptop-lipbufthick
+      if ((lipbufthick*2.0d0).gt.lipthick) &
+       write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
+      endif !lipthick.gt.0
+      write(iout,*) "bordliptop=",bordliptop
+      write(iout,*) "bordlipbot=",bordlipbot
+      write(iout,*) "bufliptop=",bufliptop
+      write(iout,*) "buflipbot=",buflipbot
+      write (iout,*) "SHIELD MODE",shield_mode
 
 !C-------------------------
       minim=(index(controlcard,'MINIMIZE').gt.0)
       if(me.eq.king.or..not.out1file) &
        write (iout,'(2a)') diagmeth(kdiag),&
         ' routine used to diagonalize matrices.'
+      if (shield_mode.gt.0) then
+      pi=3.141592d0
+!C VSolvSphere the volume of solving sphere
+!C      print *,pi,"pi"
+!C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
+!C there will be no distinction between proline peptide group and normal peptide
+!C group in case of shielding parameters
+      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
+      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
+      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
+      buff_shield=1.0d0
+      endif
       return
       end subroutine read_control
 !-----------------------------------------------------------------------------
          " stochastic forces of fully exposed sites"
          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
          do i=1,ntyp
-          write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),&
+          write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i),&
            gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
          enddo
         endif
 !      print *,"Processor",myrank," opened file IELEP" 
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,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')
+
 !      print *,"Processor",myrank," opened file ISIDEP" 
 !      print *,"Processor",myrank," opened parameter files" 
 #elif (defined G77)
       open (ielep,file=elename,status='old')
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old')
+      call getenv_loc('LIPTRANPAR',liptranname)
+      open (iliptranpar,file=liptranname,status='old')
+      call getenv_loc('TUBEPAR',tubename)
+      open (itube,file=tubename,status='old')
 #else
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
         readonly)
       open (ielep,file=elename,status='old',readonly)
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old',readonly)
+      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')
+
 #ifndef CRYST_SC
       call getenv_loc('ROTPARPDB',rotname_pdb)
       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
 
       open(irest2,file=rest2name,status='unknown')
       read(irest2,*) totT,EK,potE,totE,t_bath
+      totTafm=totT
 !      do i=1,2*nres
 ! AL 4/17/17: Now reading d_t(0,:) too
       do i=0,2*nres