boxshift &readpdb
authorCezary Czaplewski <czarek@chem.univ.gda.pl>
Wed, 1 Apr 2020 18:03:56 +0000 (20:03 +0200)
committerCezary Czaplewski <czarek@chem.univ.gda.pl>
Wed, 1 Apr 2020 18:03:56 +0000 (20:03 +0200)
20 files changed:
source/cluster/wham/src-HCD/DIMENSIONS
source/cluster/wham/src-HCD/Makefile-MPICH-ifort-okeanos
source/cluster/wham/src-HCD/readpdb.F
source/unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos
source/unres/src-HCD-5D/boxshift.f [new file with mode: 0644]
source/unres/src-HCD-5D/elecont.f
source/unres/src-HCD-5D/energy_p_new-sep_barrier.F
source/unres/src-HCD-5D/energy_p_new_barrier.F
source/unres/src-HCD-5D/gen_rand_conf.F
source/unres/src-HCD-5D/int_from_cart.F [new file with mode: 0644]
source/unres/src-HCD-5D/minimize_p.F
source/unres/src-HCD-5D/readpdb-mult.F [new file with mode: 0644]
source/unres/src-HCD-5D/readpdb.F
source/unres/src-HCD-5D/readrtns_CSA.F
source/unres/src-HCD-5D/rmscalc.F
source/unres/src-HCD-5D/sc_move.F
source/unres/src-HCD-5D/stochfric.F
source/unres/src-HCD-5D/unres.F
source/wham/src-HCD/Makefile_MPICH_ifort-okeanos
source/wham/src-HCD/readpdb.F

index 80ac845..909354c 100644 (file)
@@ -85,3 +85,6 @@ C Maximum number of SC local term fitting function coefficiants
 C Maximum number of terms in SC bond-stretching potential
       integer maxbondterm
       parameter (maxbondterm=3)
+C Maximum number of generated conformations
+      integer mxio
+      parameter (mxio=2)
index 4f7f61f..d1e64a7 100644 (file)
@@ -1,7 +1,7 @@
 #INSTALL_DIR = /opt/cray/mpt/7.3.2/gni/mpich-intel/15.0
 FC = ftn
-OPT = -O3 -ip -mcmodel=medium -shared-intel -dynamic
-#OPT = -CB -g -mcmodel=medium -shared-intel -dynamic
+#OPT = -O3 -ip -mcmodel=medium -shared-intel -dynamic
+OPT = -CB -g -mcmodel=medium -shared-intel -dynamic
 FFLAGS =  ${OPT} -c -I. -Iinclude_unres -I$(INSTALL_DIR)/include
 LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdrf.a
 
index dc6aa0a..2e73601 100644 (file)
@@ -1,9 +1,9 @@
-      subroutine readpdb(lprint)
+      subroutine readpdb
 C Read the PDB file and convert the peptide geometry into virtual-chain 
 C geometry.
-      implicit none
+      implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
-      include 'COMMON.CONTROL'
+      include 'COMMON.FRAG'
       include 'COMMON.LOCAL'
       include 'COMMON.VAR'
       include 'COMMON.CHAIN'
@@ -11,147 +11,176 @@ C geometry.
       include 'COMMON.IOUNITS'
       include 'COMMON.GEO'
       include 'COMMON.NAMES'
-      include 'COMMON.SBRIDGE'
-      character*3 seq,atom,res
+      include 'COMMON.CONTROL'
+      integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity
+      logical lprn /.false./,fail
+      double precision e1(3),e2(3),e3(3)
+      double precision dcj,efree_temp
+      character*3 seq,res
+      character*5 atom
       character*80 card
-      double precision sccor(3,50)
-      integer i,j,iii,ibeg,ishift,ishift1,ity,ires,ires_old
-      double precision dcj
-      integer rescode,kkk,lll,icha,cou,kupa,iprzes
-      logical lprint
+      double precision sccor(3,20)
+      integer rescode
+      efree_temp=0.0d0
       ibeg=1
       ishift1=0
+      ishift=0
+c      write (2,*) "UNRES_PDB",unres_pdb
+      ires=0
+      ires_old=0
+      iii=0
+      lsecondary=.false.
+      nhfrag=0
+      nbfrag=0
       do
         read (ipdbin,'(a80)',end=10) card
-        if (card(:3).eq.'END') then
-          goto 10
-        else if (card(:3).eq.'TER') then
-C End current chain
-c          ires_old=ires+1 
-          ires_old=ires+2
-          itype(ires_old-1)=ntyp1 
-          itype(ires_old)=ntyp1
-          ibeg=2
-c          write (iout,*) "Chain ended",ires,ishift,ires_old
-          call sccenter(ires,iii,sccor)
+c        write (iout,'(a)') card
+        if (card(:5).eq.'HELIX') then
+         nhfrag=nhfrag+1
+         lsecondary=.true.
+         read(card(22:25),*) hfrag(1,nhfrag)
+         read(card(34:37),*) hfrag(2,nhfrag)
+        endif
+        if (card(:5).eq.'SHEET') then
+         nbfrag=nbfrag+1
+         lsecondary=.true.
+         read(card(24:26),*) bfrag(1,nbfrag)
+         read(card(35:37),*) bfrag(2,nbfrag)
+crc----------------------------------------
+crc  to be corrected !!!
+         bfrag(3,nbfrag)=bfrag(1,nbfrag)
+         bfrag(4,nbfrag)=bfrag(2,nbfrag)
+crc----------------------------------------
         endif
+        if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
+c Read free energy
+        if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
 C Fish out the ATOM cards.
         if (index(card(1:4),'ATOM').gt.0) then  
-          read (card(14:16),'(a3)') atom
-          if (atom.eq.'CA' .or. atom.eq.'CH3') then
+          read (card(12:16),*) atom
+c          write (iout,*) "! ",atom," !",ires
+c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
+          read (card(23:26),*) ires
+          read (card(18:20),'(a3)') res
+c          write (iout,*) "ires",ires,ires-ishift+ishift1,
+c     &      " ires_old",ires_old
+c          write (iout,*) "ishift",ishift," ishift1",ishift1
+c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
+          if (ires-ishift+ishift1.ne.ires_old) then
 C Calculate the CM of the preceding residue.
+c            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
             if (ibeg.eq.0) then
-              call sccenter(ires,iii,sccor)
+c              write (iout,*) "Calculating sidechain center iii",iii
+              if (unres_pdb) then
+                do j=1,3
+                  dc(j,ires)=sccor(j,iii)
+                enddo
+              else
+                call sccenter(ires_old,iii,sccor)
+              endif
+              iii=0
             endif
 C Start new residue.
-c            write (iout,'(a80)') card
-            read (card(23:26),*) ires
-            read (card(18:20),'(a3)') res
-            if (ibeg.eq.1) then
+            if (res.eq.'Cl-' .or. res.eq.'Na+') then
+              ires=ires_old
+              cycle
+            else if (ibeg.eq.1) then
+c              write (iout,*) "BEG ires",ires
               ishift=ires-1
               if (res.ne.'GLY' .and. res.ne. 'ACE') then
                 ishift=ishift-1
                 itype(1)=ntyp1
               endif
-c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
+              ires=ires-ishift+ishift1
+              ires_old=ires
+c              write (iout,*) "ishift",ishift," ires",ires,
+c     &         " ires_old",ires_old
               ibeg=0          
-            else if (ibeg.eq.2) then
-c Start a new chain
-              ishift=-ires_old+ires-1
-c              write (iout,*) "New chain started",ires,ishift
-              ibeg=0
+            else
+              ishift=ishift-(ires-ishift+ishift1-ires_old-1)
+              ires=ires-ishift+ishift1
+              ires_old=ires
             endif
-            ires=ires-ishift
-c            write (2,*) "ires",ires," ishift",ishift
-            if (res.eq.'ACE') then
-              ity=10
+            if (res.eq.'ACE' .or. res.eq.'NHE') then
+              itype(ires)=10
             else
               itype(ires)=rescode(ires,res,0)
             endif
+          else
+            ires=ires-ishift+ishift1
+          endif
+c          write (iout,*) "ires_old",ires_old," ires",ires
+          if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+c            ishift1=ishift1+1
+          endif
+c          write (2,*) "ires",ires," res ",res," ity",ity
+          if (atom.eq.'CA' .or. atom.eq.'CH3' .or. 
+     &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
-            read(card(61:66),*) bfac(ires)
-c            write (iout,'(2i3,2x,a,3f8.3,5x,f8.3)') 
-c     &       ires,itype(ires),res,(c(j,ires),j=1,3),bfac(ires)
-            iii=1
+c            write (iout,*) "backbone ",atom 
+#ifdef DEBUG
+            write (iout,'(2i3,2x,a,3f8.3)') 
+     &      ires,itype(ires),res,(c(j,ires),j=1,3)
+#endif
+            iii=iii+1
             do j=1,3
               sccor(j,iii)=c(j,ires)
             enddo
-          else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and.
-     &             atom(1:1).ne.'Q' .and. atom(1:2).ne.'1H' .and.
-     &             atom(1:2).ne.'2H' .and. atom(1:2).ne.'3H' .and.
-     &             atom.ne.'N  ' .and. atom.ne.'C   ' .and.
-     &             atom.ne.'OXT' ) then
+            if (ishift.ne.0) then
+              ires_ca=ires+ishift-ishift1
+            else
+              ires_ca=ires
+            endif
+c            write (*,*) card(23:27),ires,itype(ires)
+          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.
+     &             atom.ne.'OXT' .and. atom(:2).ne.'3H') then
+c            write (iout,*) "sidechain ",atom
             iii=iii+1
-c            write (iout,*) res,ires,iii,atom
             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
-c            write (iout,'(3f8.3)') (sccor(j,iii),j=1,3)
           endif
         endif
       enddo
-   10 write (iout,'(a,i5)') ' Nres: ',ires
-C Calculate dummy residue coordinates inside the "chain" of a multichain
-C system
-      nres=ires
-      do i=2,nres-1
-c        write (iout,*) i,itype(i)
-
-        if (itype(i).eq.ntyp1) then
-         if (itype(i+1).eq.ntyp1) then
-C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
-C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
-C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
-C           if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the last dummy residue
-C            call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
-C            if (fail) then
-C              e2(1)=0.0d0
-C              e2(2)=1.0d0
-C              e2(3)=0.0d0
-C            endif !fail
-C            do j=1,3
-C             c(j,i)=c(j,i-1)-1.9d0*e2(j)
-C            enddo
-C           else   !unres_pdb
-           do j=1,3
-             dcj=(c(j,i-2)-c(j,i-3))/2.0
-             c(j,i)=c(j,i-1)+dcj
-             c(j,nres+i)=c(j,i)
-           enddo     
-C          endif   !unres_pdb
-         else     !itype(i+1).eq.ntyp1
-C          if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the first dummy residue
-C            call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
-C            if (fail) then
-C              e2(1)=0.0d0
-C              e2(2)=1.0d0
-C              e2(3)=0.0d0
-C            endif
-C            do j=1,3
-C              c(j,i)=c(j,i+1)-1.9d0*e2(j)
-C            enddo
-C          else !unres_pdb
-           do j=1,3
-            dcj=(c(j,i+3)-c(j,i+2))/2.0
-            c(j,i)=c(j,i+1)-dcj
-            c(j,nres+i)=c(j,i)
-           enddo
-C          endif !unres_pdb
-         endif !itype(i+1).eq.ntyp1
-        endif  !itype.eq.ntyp1
-      enddo
+   10 continue
+#ifdef DEBUG
+      write (iout,'(a,i5)') ' Number of residues found: ',ires
+#endif
+      if (ires.eq.0) return
 C Calculate the CM of the last side chain.
-      call sccenter(ires,iii,sccor)
+      if (iii.gt.0)  then
+      if (unres_pdb) then
+        do j=1,3
+          dc(j,ires)=sccor(j,iii)
+        enddo
+      else
+        call sccenter(ires,iii,sccor)
+      endif
+      endif
+      nres=ires
       nsup=nres
       nstart_sup=1
       if (itype(nres).ne.10) then
         nres=nres+1
         itype(nres)=ntyp1
+        if (unres_pdb) then
+C 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)-3.8d0*e2(j)
+          enddo
+        else
         do j=1,3
-          dcj=(c(j,nres-2)-c(j,nres-3))/2.0
+          dcj=c(j,nres-2)-c(j,nres-3)
           c(j,nres)=c(j,nres-1)+dcj
           c(j,2*nres)=c(j,nres)
         enddo
+        endif
       endif
       do i=2,nres-1
         do j=1,3
@@ -165,23 +194,57 @@ C Calculate the CM of the last side chain.
       if (itype(1).eq.ntyp1) then
         nsup=nsup-1
         nstart_sup=2
+        if (unres_pdb) then
+C 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)-3.8d0*e2(j)
+          enddo
+        else
         do j=1,3
-          dcj=(c(j,4)-c(j,3))/2.0
+          dcj=c(j,4)-c(j,3)
           c(j,1)=c(j,2)-dcj
           c(j,nres+1)=c(j,1)
         enddo
+        endif
       endif
+C Copy the coordinates to reference coordinates
+c      do i=1,2*nres
+c        do j=1,3
+c          cref(j,i)=c(j,i)
+c        enddo
+c      enddo
 C Calculate internal coordinates.
-      if (lprint) then
-      write (iout,100)
+      if (lprn) then
+      write (iout,'(/a)') 
+     &  "Cartesian coordinates of the reference structure"
+      write (iout,'(a,3(3x,a5),5x,3(3x,a5))') 
+     & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
       do ires=1,nres
+        write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') 
+     &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
+     &    (c(j,ires+nres),j=1,3)
+      enddo
+      endif
+C Calculate internal coordinates.
+       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)') 
      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
      &    (c(j,nres+ires),j=1,3)
-      enddo
-      endif
+       enddo
       call int_from_cart(.true.,.false.)
-      call flush(iout)
+      call sc_loc_geom(.false.)
+      do i=1,nres
+        thetaref(i)=theta(i)
+        phiref(i)=phi(i)
+      enddo
       do i=1,nres-1
         do j=1,3
           dc(j,i)=c(j,i+1)-c(j,i)
@@ -198,24 +261,30 @@ c     &   vbld_inv(i+nres)
       enddo
 c      call chainbuild
 C Copy the coordinates to reference coordinates
-      do i=1,nres
+      do i=1,2*nres
         do j=1,3
           cref(j,i)=c(j,i)
-          cref(j,i+nres)=c(j,i+nres)
         enddo
       enddo
-  100 format ('Residue    alpha-carbon coordinates    ',
-     &          '     centroid coordinates'/
-     1          '         ', 6X,'X',7X,'Y',7X,'Z',
-     &                          12X,'X',7X,'Y',7X,'Z')
-  110 format (a,'(',i3,')',6f12.5)
 
+
+      do j=1,nbfrag     
+        do i=1,4                                                       
+         bfrag(i,j)=bfrag(i,j)-ishift
+        enddo
+      enddo
+
+      do j=1,nhfrag
+        do i=1,2
+         hfrag(i,j)=hfrag(i,j)-ishift
+        enddo
+      enddo
       ishift_pdb=ishift
       return
       end
 c---------------------------------------------------------------------------
       subroutine int_from_cart(lside,lprn)
-      implicit none
+      implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'COMMON.LOCAL'
       include 'COMMON.VAR'
@@ -224,56 +293,61 @@ c---------------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.GEO'
       include 'COMMON.NAMES'
-      character*3 seq,atom,res
+      character*3 seq,res
+c      character*5 atom
       character*80 card
-      double precision sccor(3,50)
+      dimension sccor(3,20)
       integer rescode
-      double precision dist,alpha,beta,di
-      integer i,j,iti
       logical lside,lprn
-      if (lprn) then 
+       if (lprn) then 
         write (iout,'(/a)') 
      &  'Internal coordinates calculated from crystal structure.'
         if (lside) then 
           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
-     & '       Phi','    Dsc_id','       Dsc','     Alpha',
-     & '     Omega'
+     & '     Gamma','    Dsc_id','       Dsc','     Alpha',
+     & '     Beta '
         else 
           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
-     & '       Phi'
+     & '     Gamma'
         endif
-      endif
-      do i=2,nres
+       endif
+      do i=1,nres-1
         iti=itype(i)
-c        write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1)
-        if (itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1 .and.
-     &    (dist(i,i-1).lt.1.0D0 .or. dist(i,i-1).gt.6.0D0)) then
+        if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
-          stop
+ctest          stop
         endif
-        vbld(i)=dist(i-1,i)
-        vbld_inv(i)=1.0d0/vbld(i)
-        theta(i+1)=alpha(i-1,i,i+1)
+        vbld(i+1)=dist(i,i+1)
+        vbld_inv(i+1)=1.0d0/vbld(i+1)
+        if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
       enddo
-c      if (itype(1).eq.ntyp1) then
-c        do j=1,3
-c          c(j,1)=c(j,2)+(c(j,3)-c(j,4))
-c        enddo
-c      endif
-c      if (itype(nres).eq.ntyp1) then
-c        do j=1,3
-c          c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3))
-c        enddo
+c      if (unres_pdb) then
+c        if (itype(1).eq.ntyp1) then
+c          theta(3)=90.0d0*deg2rad
+c          phi(4)=180.0d0*deg2rad
+c          vbld(2)=3.8d0
+c          vbld_inv(2)=1.0d0/vbld(2)
+c        endif
+c        if (itype(nres).eq.ntyp1) then
+c          theta(nres)=90.0d0*deg2rad
+c          phi(nres)=180.0d0*deg2rad
+c          vbld(nres)=3.8d0
+c          vbld_inv(nres)=1.0d0/vbld(2)
+c        endif
 c      endif
       if (lside) then
         do i=2,nres-1
           do j=1,3
-            c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
+            c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
+     &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
           enddo
           iti=itype(i)
           di=dist(i,nres+i)
-           vbld(i+nres)=di
+C 10/03/12 Adam: Correction for zero SC-SC bond length
+          if (itype(i).ne.10 .and. itype(i).ne.ntyp1. and. di.eq.0.0d0)
+     &     di=dsc(itype(i))
+          vbld(i+nres)=di
           if (itype(i).ne.10) then
             vbld_inv(i+nres)=1.0d0/di
           else
@@ -283,41 +357,21 @@ c      endif
             alph(i)=alpha(nres+i,i,maxres2)
             omeg(i)=beta(nres+i,i,maxres2,i+1)
           endif
-          if (iti.ne.10) then
-            alph(i)=alpha(nres+i,i,maxres2)
-            omeg(i)=beta(nres+i,i,maxres2,i+1)
-          endif
-          if (lprn)
-     &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
-     &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,
-     &    rad2deg*alph(i),rad2deg*omeg(i)
+           if (lprn)
+     &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
+     &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
+     &     rad2deg*alph(i),rad2deg*omeg(i)
         enddo
       else if (lprn) then
         do i=2,nres
           iti=itype(i)
           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
-     &    rad2deg*theta(i),rad2deg*phi(i)
+     &     rad2deg*theta(i),rad2deg*phi(i)
         enddo
       endif
       return
       end
-c---------------------------------------------------------------------------
-      subroutine sccenter(ires,nscat,sccor)
-      implicit none
-      include 'DIMENSIONS'
-      include 'COMMON.CHAIN'
-      integer ires,nscat,i,j
-      double precision sccor(3,50),sccmj
-      do j=1,3
-        sccmj=0.0D0
-        do i=1,nscat
-          sccmj=sccmj+sccor(j,i) 
-        enddo
-        dc(j,ires)=sccmj/nscat
-      enddo
-      return
-      end
-c---------------------------------------------------------------------------
+c-------------------------------------------------------------------------------
       subroutine sc_loc_geom(lprn)
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
@@ -329,7 +383,6 @@ c---------------------------------------------------------------------------
       include 'COMMON.GEO'
       include 'COMMON.NAMES'
       include 'COMMON.CONTROL'
-      include 'COMMON.SETUP'
       double precision x_prime(3),y_prime(3),z_prime(3)
       logical lprn
       do i=1,nres-1
@@ -338,7 +391,7 @@ c---------------------------------------------------------------------------
         enddo
       enddo
       do i=2,nres-1
-        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+        if (itype(i).ne.10) then
           do j=1,3
             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
           enddo
@@ -358,7 +411,7 @@ c---------------------------------------------------------------------------
         sinfac2=0.5d0/(1.0d0-costtab(i+1))
         sinfac=dsqrt(sinfac2)
         it=itype(i)
-        if (it.ne.10 .and. itype(i).ne.ntyp1) then
+        if (it.ne.10) then
 c
 C  Compute the axes of tghe local cartesian coordinates system; store in
 c   x_prime, y_prime and z_prime 
@@ -372,10 +425,7 @@ c
           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
         enddo
-c        write (iout,*) "x_prime",(x_prime(j),j=1,3)
-c        write (iout,*) "y_prime",(y_prime(j),j=1,3)
         call vecpr(x_prime,y_prime,z_prime)
-c        write (iout,*) "z_prime",(z_prime(j),j=1,3)
 c
 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 C to local coordinate system. Store in xx, yy, zz.
@@ -399,16 +449,30 @@ c
         endif
       enddo
       if (lprn) then
-        write (iout,*) "xxref,yyref,zzref"
         do i=2,nres
           iti=itype(i)
-          write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),yyref(i),
-     &     zzref(i)
+          write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
+     &      yyref(i),zzref(i)
         enddo
       endif
       return
       end
 c---------------------------------------------------------------------------
+      subroutine sccenter(ires,nscat,sccor)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      dimension sccor(3,20)
+      do j=1,3
+        sccmj=0.0D0
+        do i=1,nscat
+          sccmj=sccmj+sccor(j,i) 
+        enddo
+        dc(j,ires)=sccmj/nscat
+      enddo
+      return
+      end
+c---------------------------------------------------------------------------
       subroutine bond_regular
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'   
@@ -420,15 +484,14 @@ c---------------------------------------------------------------------------
       do i=1,nres-1
        vbld(i+1)=vbl
        vbld_inv(i+1)=1.0d0/vbld(i+1)
-       vbld(i+1+nres)=dsc(iabs(itype(i+1)))
-       vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1)))
+       vbld(i+1+nres)=dsc(itype(i+1))
+       vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
 c       print *,vbld(i+1),vbld(i+1+nres)
       enddo
       return
       end
-c---------------------------------------------------------------------------
       subroutine readpdb_template(k)
-C Read the PDB file for read_constr_homology with read2sigma
+C Read the PDB file with gaps for read_constr_homology with read2sigma
 C and convert the peptide geometry into virtual-chain geometry.
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
@@ -440,7 +503,6 @@ C and convert the peptide geometry into virtual-chain geometry.
       include 'COMMON.GEO'
       include 'COMMON.NAMES'
       include 'COMMON.CONTROL'
-      include 'COMMON.SETUP'
       integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity
       logical lprn /.false./,fail
       double precision e1(3),e2(3),e3(3)
@@ -449,10 +511,8 @@ C and convert the peptide geometry into virtual-chain geometry.
       character*5 atom
       character*80 card
       double precision sccor(3,20)
-      integer rescode,iterter(maxres)
-      do i=1,maxres
-         iterter(i)=0
-      enddo
+      integer rescode
+      efree_temp=0.0d0
       ibeg=1
       ishift1=0
       ishift=0
@@ -463,27 +523,10 @@ c      write (2,*) "UNRES_PDB",unres_pdb
       lsecondary=.false.
       nhfrag=0
       nbfrag=0
-      do
+      do 
         read (ipdbin,'(a80)',end=10) card
-        if (card(:3).eq.'END') then
-          goto 10
-        else if (card(:3).eq.'TER') then
-C End current chain
-          ires_old=ires+2
-          itype(ires_old-1)=ntyp1 
-          iterter(ires_old-1)=1
-          itype(ires_old)=ntyp1
-          iterter(ires_old)=1
-          ibeg=2
-c          write (iout,*) "Chain ended",ires,ishift,ires_old
-          if (unres_pdb) then
-            do j=1,3
-              dc(j,ires)=sccor(j,iii)
-            enddo
-          else 
-            call sccenter(ires,iii,sccor)
-          endif
-        endif
+c        write (iout,'(a)') card
+        if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
 C Fish out the ATOM cards.
         if (index(card(1:4),'ATOM').gt.0) then  
           read (card(12:16),*) atom
@@ -497,7 +540,9 @@ c          write (iout,*) "ishift",ishift," ishift1",ishift1
 c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
           if (ires-ishift+ishift1.ne.ires_old) then
 C Calculate the CM of the preceding residue.
+c            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
             if (ibeg.eq.0) then
+c              write (iout,*) "Calculating sidechain center iii",iii
               if (unres_pdb) then
                 do j=1,3
                   dc(j,ires)=sccor(j,iii)
@@ -522,13 +567,6 @@ c              write (iout,*) "BEG ires",ires
               ires_old=ires
 c              write (iout,*) "ishift",ishift," ires",ires,
 c     &         " ires_old",ires_old
-c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
-              ibeg=0          
-            else if (ibeg.eq.2) then
-c Start a new chain
-              ishift=-ires_old+ires-1
-              ires=ires_old+1
-c              write (iout,*) "New chain started",ires,ishift
               ibeg=0          
             else
               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
@@ -544,14 +582,14 @@ c              write (iout,*) "New chain started",ires,ishift
             ires=ires-ishift+ishift1
           endif
 c          write (iout,*) "ires_old",ires_old," ires",ires
-c          if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+          if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 c            ishift1=ishift1+1
-c          endif
+          endif
 c          write (2,*) "ires",ires," res ",res," ity",ity
           if (atom.eq.'CA' .or. atom.eq.'CH3' .or. 
      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
-c            write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
+c            write (iout,*) "backbone ",atom 
 #ifdef DEBUG
             write (iout,'(2i3,2x,a,3f8.3)') 
      &      ires,itype(ires),res,(c(j,ires),j=1,3)
@@ -576,60 +614,13 @@ c            write (iout,*) "sidechain ",atom
           endif
         endif
       enddo
-   10 write (iout,'(a,i5)') ' Nres: ',ires
-C Calculate dummy residue coordinates inside the "chain" of a multichain
-C system
-      nres=ires
-      do i=2,nres-1
-c        write (iout,*) i,itype(i),itype(i+1)
-        if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
-         if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
-C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
-C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
-C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
-           if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the last dummy residue
-            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
-            do j=1,3
-             c(j,i)=c(j,i-1)-1.9d0*e2(j)
-            enddo
-           else   !unres_pdb
-           do j=1,3
-             dcj=(c(j,i-2)-c(j,i-3))/2.0
-            if (dcj.eq.0) dcj=1.23591524223
-             c(j,i)=c(j,i-1)+dcj
-             c(j,nres+i)=c(j,i)
-           enddo     
-          endif   !unres_pdb
-         else     !itype(i+1).eq.ntyp1
-          if (unres_pdb) then
-C 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
-            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+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 !itype(i+1).eq.ntyp1
-        endif  !itype.eq.ntyp1
-      enddo
+   10 continue
+#ifdef DEBUG
+      write (iout,'(a,i5)') ' Number of residues found: ',ires
+#endif
+      if (ires.eq.0) return
 C Calculate the CM of the last side chain.
+      if (iii.gt.0)  then
       if (unres_pdb) then
         do j=1,3
           dc(j,ires)=sccor(j,iii)
@@ -637,31 +628,19 @@ C Calculate the CM of the last side chain.
       else
         call sccenter(ires,iii,sccor)
       endif
+      endif
+      nres=ires
       nsup=nres
       nstart_sup=1
       if (itype(nres).ne.10) then
         nres=nres+1
         itype(nres)=ntyp1
-        if (unres_pdb) then
-C 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
         do j=1,3
-          dcj=(c(j,nres-2)-c(j,nres-3))/2.0
-            if (dcj.eq.0) dcj=1.23591524223
+          dcj=c(j,nres-2)-c(j,nres-3)
           c(j,nres)=c(j,nres-1)+dcj
           c(j,2*nres)=c(j,nres)
         enddo
       endif
-      endif
       do i=2,nres-1
         do j=1,3
           c(j,i+nres)=dc(j,i)
@@ -683,24 +662,18 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue
             e2(3)=0.0d0
           endif
           do j=1,3
-            c(j,1)=c(j,2)-1.9d0*e2(j)
+            c(j,1)=c(j,2)-3.8d0*e2(j)
           enddo
         else
         do j=1,3
-          dcj=(c(j,4)-c(j,3))/2.0
+          dcj=c(j,4)-c(j,3)
           c(j,1)=c(j,2)-dcj
           c(j,nres+1)=c(j,1)
         enddo
         endif
       endif
-C Copy the coordinates to reference coordinates
-c      do i=1,2*nres
-c        do j=1,3
-c          cref(j,i)=c(j,i)
-c        enddo
-c      enddo
 C Calculate internal coordinates.
-      if (out_template_coord) then
+      if (lprn) then
       write (iout,'(/a)') 
      &  "Cartesian coordinates of the reference structure"
       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') 
@@ -712,9 +685,15 @@ C Calculate internal coordinates.
       enddo
       endif
 C Calculate internal coordinates.
-c      call int_from_cart1(.false.)
-      call int_from_cart(.true.,.true.)
-      call sc_loc_geom(.true.)
+       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)') 
+     &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
+     &    (c(j,nres+ires),j=1,3)
+       enddo
+      call int_from_cart(.true.,.false.)
+      call sc_loc_geom(.false.)
       do i=1,nres
         thetaref(i)=theta(i)
         phiref(i)=phi(i)
@@ -733,19 +712,16 @@ c      call int_from_cart1(.false.)
 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 c     &   vbld_inv(i+nres)
       enddo
-      do i=1,nres
-        do j=1,3
-          cref(j,i)=c(j,i)
-          cref(j,i+nres)=c(j,i+nres)
-        enddo
-      enddo
+c      call chainbuild
+C Copy the coordinates to reference coordinates
       do i=1,2*nres
         do j=1,3
+          cref(j,i)=c(j,i)
           chomo(j,i,k)=c(j,i)
         enddo
       enddo
 
+
+      ishift_pdb=ishift
       return
       end
-      
-
index b7c9d19..9826556 100644 (file)
@@ -40,7 +40,8 @@ object = unres.o arcos.o cartprint.o chainbuild.o convert.o initialize_p.o \
         cart2intgrad.o checkder_p.o contact_cp econstr_local.o econstr_qlike.o \
        econstrq-PMF.o PMFprocess.o energy_p_new_barrier.o make_xx_list.o \
        energy_p_new-sep_barrier.o gradient_p.o minimize_p.o sumsld.o \
-        cored.o rmdd.o geomout.o readpdb.o regularize.o thread.o fitsq.o mcm.o \
+        cored.o rmdd.o geomout.o readpdb.o int_from_cart.o regularize.o \
+       thread.o fitsq.o mcm.o \
         mc.o bond_move.o refsys.o check_sc_distr.o check_bond.o contact.o \
         eigen.o blas.o add.o entmcm.o minim_mcmf.o \
         together.o csa.o minim_jlee.o shift.o diff12.o bank.o newconf.o ran.o \
@@ -84,6 +85,15 @@ E0LL2Y: ${object} xdrf/libxdrf.a
        ${FC} ${FFLAGS} cinfo.f
        ${FC} ${OPT} ${object} cinfo.o ${LIBS}  -o ${BIN}
 
+E0LL2Y_DFA: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI \
+       -DSPLITELE -DLANG0 -DFOURBODY -DDFA
+E0LL2Y_DFA: BIN = ~/bin/unres_ifort_MPICH-okeanos_E0LL2Y-HCD-DFA.exe
+E0LL2Y_DFA: ${object} dfa.o xdrf/libxdrf.a
+       gcc -o compinfo compinfo.c
+       ./compinfo | true
+       ${FC} ${FFLAGS} cinfo.f
+       ${FC} ${OPT} ${object} dfa.o cinfo.o ${LIBS}  -o ${BIN}
+
 NEWCORR: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI \
        -DSPLITELE -DLANG0 -DNEWCORR -DCORRCD #-DFOURBODY #-DMYGAUSS #-DTIMING
 NEWCORR: BIN = ~/bin/unres_ifort_MPICH-okeanos_SC-HCD.exe
@@ -151,6 +161,9 @@ cartder.o : cartder.F
 readpdb.o : readpdb.F
        ${FC} ${FFLAGS2} ${CPPFLAGS} readpdb.F
 
+readpdb.o : readpdb.F
+       ${FC} ${FFLAGS2} ${CPPFLAGS} readpdb-mult.F
+
 sumsld.o : sumsld.f
        ${FC} ${FFLAGS} ${CPPFLAGS} sumsld.f
         
diff --git a/source/unres/src-HCD-5D/boxshift.f b/source/unres/src-HCD-5D/boxshift.f
new file mode 100644 (file)
index 0000000..29d3406
--- /dev/null
@@ -0,0 +1,101 @@
+
+c------------------------------------------------------------------------
+      double precision function boxshift(x,boxsize)
+      implicit none
+      double precision x,boxsize
+      double precision xtemp
+      xtemp=dmod(x,boxsize)
+      if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then
+        boxshift=xtemp-boxsize
+      else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then
+        boxshift=xtemp+boxsize
+      else
+        boxshift=xtemp
+      endif
+      return
+      end
+c--------------------------------------------------------------------------
+      subroutine closest_img(xi,yi,zi,xj,yj,zj)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      integer xshift,yshift,zshift,subchap
+      double precision dist_init,xj_safe,yj_safe,zj_safe,
+     & xj_temp,yj_temp,zj_temp,dist_temp
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+      subchap=0
+      do xshift=-1,1
+        do yshift=-1,1
+          do zshift=-1,1
+            xj=xj_safe+xshift*boxxsize
+            yj=yj_safe+yshift*boxysize
+            zj=zj_safe+zshift*boxzsize
+            dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+            if(dist_temp.lt.dist_init) then
+              dist_init=dist_temp
+              xj_temp=xj
+              yj_temp=yj
+              zj_temp=zj
+              subchap=1
+            endif
+          enddo
+        enddo
+      enddo
+      if (subchap.eq.1) then
+        xj=xj_temp-xi
+        yj=yj_temp-yi
+        zj=zj_temp-zi
+      else
+        xj=xj_safe-xi
+        yj=yj_safe-yi
+        zj=zj_safe-zi
+      endif
+      return
+      end
+c--------------------------------------------------------------------------
+      subroutine to_box(xi,yi,zi)
+      implicit none
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      double precision xi,yi,zi
+      xi=dmod(xi,boxxsize)
+      if (xi.lt.0.0d0) xi=xi+boxxsize
+      yi=dmod(yi,boxysize)
+      if (yi.lt.0.0d0) yi=yi+boxysize
+      zi=dmod(zi,boxzsize)
+      if (zi.lt.0.0d0) zi=zi+boxzsize
+      return
+      end
+c--------------------------------------------------------------------------
+      subroutine lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+      implicit none
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      double precision xi,yi,zi,sslipi,ssgradlipi
+      double precision fracinbuf
+      double precision sscalelip,sscagradlip
+
+      if ((zi.gt.bordlipbot).and.(zi.lt.bordliptop)) then
+C the energy transfer exist
+        if (zi.lt.buflipbot) then
+C what fraction I am in
+          fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+          sslipi=sscalelip(fracinbuf)
+          ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+        elseif (zi.gt.bufliptop) then
+          fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
+          sslipi=sscalelip(fracinbuf)
+          ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+        else
+          sslipi=1.0d0
+          ssgradlipi=0.0
+        endif
+      else
+        sslipi=0.0d0
+        ssgradlipi=0.0
+      endif
+      return
+      end
index 690fd44..bf9056a 100644 (file)
@@ -12,7 +12,7 @@
       double precision app_(2,2),bpp_(2,2),rpp_(2,2)
       integer ncont,icont(2,maxcont)
       double precision econt(maxcont)
-      integer xshift,yshift,zshift
+      double precision boxshift
 *
 * Load the constants of peptide bond - peptide bond interactions.
 * Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
@@ -53,13 +53,8 @@ c      data epp    / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
         xmedi=xi+0.5*dxi
         ymedi=yi+0.5*dyi
         zmedi=zi+0.5*dzi
+        call to_box(xmedi,ymedi,zmedi)
 c        write (iout,*) "i",xmedi,ymedi,zmedi
-          xmedi=mod(xmedi,boxxsize)
-          if (xmedi.lt.0) xmedi=xmedi+boxxsize
-          ymedi=mod(ymedi,boxysize)
-          if (ymedi.lt.0) ymedi=ymedi+boxysize
-          zmedi=mod(zmedi,boxzsize)
-          if (zmedi.lt.0) zmedi=zmedi+boxzsize
 c        write (iout,*) "i",xmedi,ymedi,zmedi
         do 4 j=i+2,nct-1
 c          write (iout,*) "i",i," j",j
@@ -80,47 +75,10 @@ c          write (iout,*) "i",i," j",j
           yj=c(2,j)+0.5*dyj
           zj=c(3,j)+0.5*dzj
 c          write (iout,*) "j",xj,yj,zj
-          xj=mod(xj,boxxsize)
-          if (xj.lt.0) xj=xj+boxxsize
-          yj=mod(yj,boxysize)
-          if (yj.lt.0) yj=yj+boxysize
-          zj=mod(zj,boxzsize)
-          if (zj.lt.0) zj=zj+boxzsize
-c          write (iout,*) "j",xj,yj,zj
-      dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
-c      write (iout,*) "dist",dsqrt(dist_init)
-      xj_safe=xj
-      yj_safe=yj
-      zj_safe=zj
-      isubchap=0
-      do xshift=-1,1
-      do yshift=-1,1
-      do zshift=-1,1
-          xj=xj_safe+xshift*boxxsize
-          yj=yj_safe+yshift*boxysize
-          zj=zj_safe+zshift*boxzsize
-          dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
-c          write (iout,*) "shift",xshift,yshift,zshift," dist_temp",
-c     &      dist_temp," dist_init",dist_init
-          if(dist_temp.lt.dist_init) then
-            dist_init=dist_temp
-            xj_temp=xj
-            yj_temp=yj
-            zj_temp=zj
-            isubchap=1
-          endif
-       enddo
-       enddo
-       enddo
-       if (isubchap.eq.1) then
-          xj=xj_temp-xmedi
-          yj=yj_temp-ymedi
-          zj=zj_temp-zmedi
-       else
-          xj=xj_safe-xmedi
-          yj=yj_safe-ymedi
-          zj=zj_safe-zmedi
-       endif
+          call to_box(xj,yj,zj)
+          xj=boxshift(xj-xi,boxxsize)
+          yj=boxshift(yj-yi,boxysize)
+          zj=boxshift(zj-zi,boxzsize)
           rij=xj*xj+yj*yj+zj*zj
           rrmij=1.0/(xj*xj+yj*yj+zj*zj)
           rmij=sqrt(rrmij)
index c6c6832..96f7777 100644 (file)
@@ -85,6 +85,7 @@ c      include 'COMMON.CONTACTS'
       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
      & sigij,r0ij,rcut,sss1,sssgrad1,sqrij
       double precision sscale,sscagrad
+      double precision boxshift
 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
 c      do i=iatsc_s,iatsc_e
@@ -97,6 +98,7 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
 C
 C Calculate SC interaction energy.
 C
@@ -106,9 +108,13 @@ cd   &                  'iend=',iend(i,iint)
 c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             rij=xj*xj+yj*yj+zj*zj
             sqrij=dsqrt(rrij)
             eps0ij=eps(itypi,itypj)
@@ -187,6 +193,7 @@ c      include 'COMMON.CONTACTS'
       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
      & sigij,r0ij,rcut,sqrij,sss1,sssgrad1
       double precision sscale,sscagrad
+      double precision boxshift
 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
 c      do i=iatsc_s,iatsc_e
@@ -199,6 +206,7 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
 C Change 12/1/95
         num_conti=0
 C
@@ -210,9 +218,13 @@ cd   &                  'iend=',iend(i,iint)
 c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
 C Change 12/1/95 to calculate four-body interactions
             rij=xj*xj+yj*yj+zj*zj
             sqrij=dsqrt(rij)
@@ -285,6 +297,7 @@ C
      & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
       logical scheck
       double precision sscale,sscagrad
+      double precision boxshift
 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
 c      do i=iatsc_s,iatsc_e
@@ -297,6 +310,7 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
 C
 C Calculate SC interaction energy.
 C
@@ -304,9 +318,13 @@ c        do iint=1,nint_gr(i)
 c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             fac_augm=rrij**expon
             e_augm=augm(itypi,itypj)*fac_augm
@@ -384,6 +402,7 @@ C
      & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
       logical scheck
       double precision sscale,sscagrad
+      double precision boxshift
 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
 c      do i=iatsc_s,iatsc_e
@@ -396,6 +415,7 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
 C
 C Calculate SC interaction energy.
 C
@@ -403,9 +423,13 @@ c        do iint=1,nint_gr(i)
 c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             fac_augm=rrij**expon
             e_augm=augm(itypi,itypj)*fac_augm
@@ -478,6 +502,7 @@ C
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
       double precision sss1,sssgrad1
       double precision sscale,sscagrad
+      double precision boxshift
 c     double precision rrsave(maxdim)
       logical lprn
       evdw=0.0D0
@@ -499,6 +524,7 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -523,9 +549,13 @@ c            dscj_inv=dsc_inv(itypj)
             alf1=alp(itypi)
             alf2=alp(itypj)
             alf12=0.5D0*(alf1+alf2)
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
@@ -604,6 +634,7 @@ C
       integer itypi,itypj,itypi1,iint,ind,ikont
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
       double precision sscale,sscagrad
+      double precision boxshift
 c     double precision rrsave(maxdim)
       logical lprn
       evdw=0.0D0
@@ -625,6 +656,7 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -649,9 +681,13 @@ c            dscj_inv=dsc_inv(itypj)
             alf1=alp(itypi)
             alf2=alp(itypj)
             alf12=0.5D0*(alf1+alf2)
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
@@ -731,6 +767,7 @@ C
      & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
       double precision dist,sscale,sscagrad,sscagradlip,sscalelip
       double precision subchap,sss1,sssgrad1
+      double precision boxshift
       evdw=0.0D0
 ccccc      energy_dec=.false.
 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -748,12 +785,8 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
-          xi=mod(xi,boxxsize)
-          if (xi.lt.0) xi=xi+boxxsize
-          yi=mod(yi,boxysize)
-          if (yi.lt.0) yi=yi+boxysize
-          zi=mod(zi,boxzsize)
-          if (zi.lt.0) zi=zi+boxzsize
+        call to_box(xi,yi,zi)
+        call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -787,79 +820,27 @@ c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
             xj=c(1,nres+j)
             yj=c(2,nres+j)
             zj=c(3,nres+j)
-            xj=mod(xj,boxxsize)
-            if (xj.lt.0) xj=xj+boxxsize
-            yj=mod(yj,boxysize)
-            if (yj.lt.0) yj=yj+boxysize
-            zj=mod(zj,boxzsize)
-            if (zj.lt.0) zj=zj+boxzsize
-            if ((zj.gt.bordlipbot).and.(zj.lt.bordliptop)) then
-C the energy transfer exist
-              if (zj.lt.buflipbot) then
-C what fraction I am in
-                fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
-                sslipj=sscalelip(fracinbuf)
-                ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
-              else if (zi.gt.bufliptop) then
-                fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
-                sslipj=sscalelip(fracinbuf)
-                ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
-             else
-               sslipj=1.0d0
-               ssgradlipj=0.0
-             endif
-           else
-             sslipj=0.0d0
-             ssgradlipj=0.0
-           endif
-           aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
-     &       +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
-           bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
-     &       +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
-           dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-           xj_safe=xj
-           yj_safe=yj
-           zj_safe=zj
-           subchap=0
-           do xshift=-1,1
-             do yshift=-1,1
-               do zshift=-1,1
-                 xj=xj_safe+xshift*boxxsize
-                 yj=yj_safe+yshift*boxysize
-                 zj=zj_safe+zshift*boxzsize
-                 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-                 if (dist_temp.lt.dist_init) then
-                   dist_init=dist_temp
-                   xj_temp=xj
-                   yj_temp=yj
-                   zj_temp=zj
-                   subchap=1
-                 endif
-               enddo
-             enddo
-           enddo
-           if (subchap.eq.1) then
-             xj=xj_temp-xi
-             yj=yj_temp-yi
-             zj=zj_temp-zi
-           else
-             xj=xj_safe-xi
-             yj=yj_safe-yi
-             zj=zj_safe-zi
-           endif
-           dxj=dc_norm(1,nres+j)
-           dyj=dc_norm(2,nres+j)
-           dzj=dc_norm(3,nres+j)
-           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
-           rij=dsqrt(rrij)
-           sss1=sscale(1.0d0/rij,r_cut_int)
-           if (sss1.eq.0.0d0) cycle
-           sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
-           if (sss.lt.1.0d0) then
+            call to_box(xj,yj,zj)
+            call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+            aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &        +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
+            bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &        +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
+            dxj=dc_norm(1,nres+j)
+            dyj=dc_norm(2,nres+j)
+            dzj=dc_norm(3,nres+j)
+            rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+            rij=dsqrt(rrij)
+            sss1=sscale(1.0d0/rij,r_cut_int)
+            if (sss1.eq.0.0d0) cycle
+            sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
+            if (sss.lt.1.0d0) then
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
-             sssgrad=
+              sssgrad=
      &         sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
               sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
               call sc_angular
@@ -946,15 +927,13 @@ C
       include 'COMMON.CONTROL'
       include "COMMON.SPLITELE"
       logical lprn
-      integer xshift,yshift,zshift
       double precision evdw
       integer itypi,itypj,itypi1,iint,ind,ikont
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
       double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
-     & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
-     & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
+     & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip
       double precision dist,sscale,sscagrad,sscagradlip,sscalelip
-      double precision subchap
+      double precision boxshift
       evdw=0.0D0
 ccccc      energy_dec=.false.
 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -972,12 +951,8 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
-          xi=mod(xi,boxxsize)
-          if (xi.lt.0) xi=xi+boxxsize
-          yi=mod(yi,boxysize)
-          if (yi.lt.0) yi=yi+boxysize
-          zi=mod(zi,boxzsize)
-          if (zi.lt.0) zi=zi+boxzsize
+        call to_box(xi,yi,zi)
+        call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -1011,67 +986,15 @@ c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
             xj=c(1,nres+j)
             yj=c(2,nres+j)
             zj=c(3,nres+j)
-            xj=mod(xj,boxxsize)
-            if (xj.lt.0) xj=xj+boxxsize
-            yj=mod(yj,boxysize)
-            if (yj.lt.0) yj=yj+boxysize
-            zj=mod(zj,boxzsize)
-            if (zj.lt.0) zj=zj+boxzsize
-            if ((zj.gt.bordlipbot).and.(zj.lt.bordliptop)) then
-C the energy transfer exist
-              if (zj.lt.buflipbot) then
-C what fraction I am in
-                fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
-                sslipj=sscalelip(fracinbuf)
-                ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
-              elseif (zi.gt.bufliptop) then
-                fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
-                sslipj=sscalelip(fracinbuf)
-                ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
-              else
-                sslipj=1.0d0
-                ssgradlipj=0.0
-              endif
-            else
-              sslipj=0.0d0
-              ssgradlipj=0.0
-            endif
+            call to_box(xj,yj,zj)
+            call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
             aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
      &        +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
             bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
      &        +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
-            dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-            xj_safe=xj
-            yj_safe=yj
-            zj_safe=zj
-            subchap=0
-            do xshift=-1,1
-              do yshift=-1,1
-                do zshift=-1,1
-                  xj=xj_safe+xshift*boxxsize
-                  yj=yj_safe+yshift*boxysize
-                  zj=zj_safe+zshift*boxzsize
-                  dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-                  if(dist_temp.lt.dist_init) then
-                    dist_init=dist_temp
-                    xj_temp=xj
-                    yj_temp=yj
-                    zj_temp=zj
-                    subchap=1
-                  endif
-                enddo
-              enddo
-            enddo
-            if (subchap.eq.1) then
-              xj=xj_temp-xi
-              yj=yj_temp-yi
-              zj=zj_temp-zi
-            else
-              xj=xj_safe-xi
-              yj=yj_safe-yi
-              zj=zj_safe-zi
-            endif
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
@@ -1191,6 +1114,8 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
+        call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -1217,9 +1142,18 @@ c            dscj_inv=dsc_inv(itypj)
             alf1=alp(itypi)
             alf2=alp(itypj)
             alf12=0.5D0*(alf1+alf2)
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+            aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+            bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &       +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
@@ -1317,6 +1251,7 @@ C
      & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
      & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
       double precision dist,sscale,sscagrad,sscagradlip,sscalelip
+      double precision boxshift
       evdw=0.0D0
 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
@@ -1336,6 +1271,8 @@ c      do i=iatsc_s,iatsc_e
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
+        call to_box(xi,yi,zi)
+        call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
 c        dsci_inv=dsc_inv(itypi)
         dsci_inv=vbld_inv(i+nres)
 C
@@ -1359,9 +1296,18 @@ c            dscj_inv=dsc_inv(itypj)
             alf1=alp(itypi)
             alf2=alp(itypj)
             alf12=0.5D0*(alf1+alf2)
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+            aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+            bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &       +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
@@ -1612,12 +1558,7 @@ C     &  .or. itype(i+4).eq.ntyp1
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
-          xmedi=mod(xmedi,boxxsize)
-          if (xmedi.lt.0) xmedi=xmedi+boxxsize
-          ymedi=mod(ymedi,boxysize)
-          if (ymedi.lt.0) ymedi=ymedi+boxysize
-          zmedi=mod(zmedi,boxzsize)
-          if (zmedi.lt.0) zmedi=zmedi+boxzsize
+        call to_box(xmedi,ymedi,zmedi)
         num_conti=0
         call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
@@ -1641,12 +1582,7 @@ C     &    .or. itype(i-1).eq.ntyp1
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
-          xmedi=mod(xmedi,boxxsize)
-          if (xmedi.lt.0) xmedi=xmedi+boxxsize
-          ymedi=mod(ymedi,boxysize)
-          if (ymedi.lt.0) ymedi=ymedi+boxysize
-          zmedi=mod(zmedi,boxzsize)
-          if (zmedi.lt.0) zmedi=zmedi+boxzsize
+        call to_box(xmedi,ymedi,zmedi)
 #ifdef FOURBODY
         num_conti=num_cont_hb(i)
 #endif
@@ -1677,12 +1613,7 @@ C     &  .or. itype(i-1).eq.ntyp1
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
-          xmedi=mod(xmedi,boxxsize)
-          if (xmedi.lt.0) xmedi=xmedi+boxxsize
-          ymedi=mod(ymedi,boxysize)
-          if (ymedi.lt.0) ymedi=ymedi+boxysize
-          zmedi=mod(zmedi,boxzsize)
-          if (zmedi.lt.0) zmedi=zmedi+boxzsize
+        call to_box(xmedi,ymedi,zmedi)
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend
 #ifdef FOURBODY
         num_conti=num_cont_hb(i)
@@ -1764,6 +1695,7 @@ C-------------------------------------------------------------------------------
       double precision sss1,sssgrad1
       double precision sscale,sscagrad
       double precision scalar
+      double precision boxshift
 
 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
 #ifdef MOMENT
@@ -1796,44 +1728,10 @@ C      print *,"WCHODZE2"
       xj=c(1,j)+0.5D0*dxj
       yj=c(2,j)+0.5D0*dyj
       zj=c(3,j)+0.5D0*dzj
-      xj=mod(xj,boxxsize)
-      if (xj.lt.0) xj=xj+boxxsize
-      yj=mod(yj,boxysize)
-      if (yj.lt.0) yj=yj+boxysize
-      zj=mod(zj,boxzsize)
-      if (zj.lt.0) zj=zj+boxzsize
-      dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
-      xj_safe=xj
-      yj_safe=yj
-      zj_safe=zj
-      isubchap=0
-      do xshift=-1,1
-      do yshift=-1,1
-      do zshift=-1,1
-          xj=xj_safe+xshift*boxxsize
-          yj=yj_safe+yshift*boxysize
-          zj=zj_safe+zshift*boxzsize
-          dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
-          if(dist_temp.lt.dist_init) then
-            dist_init=dist_temp
-            xj_temp=xj
-            yj_temp=yj
-            zj_temp=zj
-            isubchap=1
-          endif
-      enddo
-      enddo
-      enddo
-      if (isubchap.eq.1) then
-         xj=xj_temp-xmedi
-         yj=yj_temp-ymedi
-         zj=zj_temp-zmedi
-      else
-         xj=xj_safe-xmedi
-         yj=yj_safe-ymedi
-         zj=zj_safe-zmedi
-      endif
-
+      call to_box(xj,yj,zj)
+      xj=boxshift(xj-xmedi,boxxsize)
+      yj=boxshift(yj-ymedi,boxysize)
+      zj=boxshift(zj-zmedi,boxzsize)
       rij=xj*xj+yj*yj+zj*zj
       rrmij=1.0D0/rij
       rij=dsqrt(rij)
@@ -2734,6 +2632,7 @@ c      write (iout,*) "evdwpp_short"
       double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
      & dist_temp, dist_init,sss_grad
       double precision sscale,sscagrad
+      double precision boxshift
       integer ikont
       evdw1=0.0D0
 C      print *,"WCHODZE"
@@ -2754,12 +2653,7 @@ c      do i=iatel_s_vdw,iatel_e_vdw
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
-        xmedi=mod(xmedi,boxxsize)
-        if (xmedi.lt.0.0d0) xmedi=xmedi+boxxsize
-        ymedi=mod(ymedi,boxysize)
-        if (ymedi.lt.0.0d0) ymedi=ymedi+boxysize
-        zmedi=mod(zmedi,boxzsize)
-        if (zmedi.lt.0.0d0) zmedi=zmedi+boxzsize
+        call to_box(xmedi,ymedi,zmedi)
         num_conti=0
 c        write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
 c     &   ' ielend',ielend_vdw(i)
@@ -2781,43 +2675,10 @@ c        do j=ielstart_vdw(i),ielend_vdw(i)
           xj=c(1,j)+0.5D0*dxj
           yj=c(2,j)+0.5D0*dyj
           zj=c(3,j)+0.5D0*dzj
-          xj=mod(xj,boxxsize)
-          if (xj.lt.0) xj=xj+boxxsize
-          yj=mod(yj,boxysize)
-          if (yj.lt.0) yj=yj+boxysize
-          zj=mod(zj,boxzsize)
-          if (zj.lt.0) zj=zj+boxzsize
-          dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
-          xj_safe=xj
-          yj_safe=yj
-          zj_safe=zj
-          isubchap=0
-          do xshift=-1,1
-          do yshift=-1,1
-          do zshift=-1,1
-              xj=xj_safe+xshift*boxxsize
-              yj=yj_safe+yshift*boxysize
-              zj=zj_safe+zshift*boxzsize
-              dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
-              if(dist_temp.lt.dist_init) then
-                dist_init=dist_temp
-                xj_temp=xj
-                yj_temp=yj
-                zj_temp=zj
-                isubchap=1
-              endif
-           enddo
-           enddo
-           enddo
-           if (isubchap.eq.1) then
-              xj=xj_temp-xmedi
-              yj=yj_temp-ymedi
-              zj=zj_temp-zmedi
-           else
-              xj=xj_safe-xmedi
-              yj=yj_safe-ymedi
-              zj=zj_safe-zmedi
-           endif
+          call to_box(xj,yj,zj)
+          xj=boxshift(xj-xmedi,boxxsize)
+          yj=boxshift(yj-ymedi,boxysize)
+          zj=boxshift(zj-zmedi,boxzsize)
           rij=xj*xj+yj*yj+zj*zj
           rrmij=1.0D0/rij
           rij=dsqrt(rij)
       logical lprint_short
       common /shortcheck/ lprint_short
       double precision ggg(3)
-      integer xshift,yshift,zshift
       integer i,iint,j,k,iteli,itypj,subchap
       double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
      & fac,e1,e2,rij
       double precision evdw2,evdw2_14,evdwij
-      double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
-     & dist_temp, dist_init
       double precision sscale,sscagrad
+      double precision boxshift
       integer ikont
       if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb
       evdw2=0.0D0
@@ -2907,12 +2766,7 @@ c      do i=iatscp_s,iatscp_e
         xi=0.5D0*(c(1,i)+c(1,i+1))
         yi=0.5D0*(c(2,i)+c(2,i+1))
         zi=0.5D0*(c(3,i)+c(3,i+1))
-        xi=mod(xi,boxxsize)
-        if (xi.lt.0) xi=xi+boxxsize
-        yi=mod(yi,boxysize)
-        if (yi.lt.0) yi=yi+boxysize
-        zi=mod(zi,boxzsize)
-        if (zi.lt.0) zi=zi+boxzsize
+        call to_box(xi,yi,zi)
 
 c        do iint=1,nscp_gr(i)
 
@@ -2928,44 +2782,10 @@ C Uncomment following three lines for Ca-p interactions
           yj=c(2,j)
           zj=c(3,j)
 c corrected by AL
-          xj=mod(xj,boxxsize)
-          if (xj.lt.0) xj=xj+boxxsize
-          yj=mod(yj,boxysize)
-          if (yj.lt.0) yj=yj+boxysize
-          zj=mod(zj,boxzsize)
-          if (zj.lt.0) zj=zj+boxzsize
-c end correction
-          dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-          xj_safe=xj
-          yj_safe=yj
-          zj_safe=zj
-          subchap=0
-          do xshift=-1,1
-          do yshift=-1,1
-          do zshift=-1,1
-          xj=xj_safe+xshift*boxxsize
-          yj=yj_safe+yshift*boxysize
-          zj=zj_safe+zshift*boxzsize
-          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-          if(dist_temp.lt.dist_init) then
-            dist_init=dist_temp
-            xj_temp=xj
-            yj_temp=yj
-            zj_temp=zj
-            subchap=1
-          endif
-          enddo
-          enddo
-          enddo
-          if (subchap.eq.1) then
-            xj=xj_temp-xi
-            yj=yj_temp-yi
-            zj=zj_temp-zi
-          else
-            xj=xj_safe-xi
-            yj=yj_safe-yi
-            zj=zj_safe-zi
-          endif
+          call to_box(xj,yj,zj)
+          xj=boxshift(xj-xi,boxxsize)
+          yj=boxshift(yj-yi,boxysize)
+          zj=boxshift(zj-zi,boxzsize)
 
           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 
@@ -3063,6 +2883,7 @@ C
      & dist_temp, dist_init
       double precision ggg(3)
       double precision sscale,sscagrad
+      double precision boxshift
       evdw2=0.0D0
       evdw2_14=0.0d0
 cd    print '(a)','Enter ESCP'
@@ -3076,12 +2897,7 @@ c     & ' iatscp_e=',iatscp_e
         xi=0.5D0*(c(1,i)+c(1,i+1))
         yi=0.5D0*(c(2,i)+c(2,i+1))
         zi=0.5D0*(c(3,i)+c(3,i+1))
-        xi=mod(xi,boxxsize)
-        if (xi.lt.0) xi=xi+boxxsize
-        yi=mod(yi,boxysize)
-        if (yi.lt.0) yi=yi+boxysize
-        zi=mod(zi,boxzsize)
-        if (zi.lt.0) zi=zi+boxzsize
+        call to_box(xi,yi,zi)
 
 c        if (lprint_short) 
 c     &    write (iout,*) "i",i," itype",itype(i),itype(i+1),
@@ -3102,49 +2918,13 @@ C Uncomment following three lines for Ca-p interactions
           yj=c(2,j)
           zj=c(3,j)
 c corrected by AL
-          xj=mod(xj,boxxsize)
-          if (xj.lt.0) xj=xj+boxxsize
-          yj=mod(yj,boxysize)
-          if (yj.lt.0) yj=yj+boxysize
-          zj=mod(zj,boxzsize)
-          if (zj.lt.0) zj=zj+boxzsize
-c end correction
-          dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-c          if (lprint_short) then
-c            write (iout,*) i,j,xi,yi,zi,xj,yj,zj
-c            write (iout,*) "dist_init",dsqrt(dist_init)
-c          endif
-          xj_safe=xj
-          yj_safe=yj
-          zj_safe=zj
-          subchap=0
-          do xshift=-1,1
-          do yshift=-1,1
-          do zshift=-1,1
-          xj=xj_safe+xshift*boxxsize
-          yj=yj_safe+yshift*boxysize
-          zj=zj_safe+zshift*boxzsize
-          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-          if(dist_temp.lt.dist_init) then
-            dist_init=dist_temp
-            xj_temp=xj
-            yj_temp=yj
-            zj_temp=zj
-            subchap=1
-          endif
-          enddo
-          enddo
-          enddo
-c          if (lprint_short) write (iout,*) "dist_temp",dsqrt(dist_temp)
-          if (subchap.eq.1) then
-             xj=xj_temp-xi
-             yj=yj_temp-yi
-             zj=zj_temp-zi
-          else
-             xj=xj_safe-xi
-             yj=yj_safe-yi
-             zj=zj_safe-zi
-          endif
+          xj=c(1,j)
+          yj=c(2,j)
+          zj=c(3,j)
+          call to_box(xj,yj,zj)
+          xj=boxshift(xj-xi,boxxsize)
+          yj=boxshift(yj-yi,boxysize)
+          zj=boxshift(zj-zi,boxzsize)
           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 c          sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
 c          sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
index ba7cbd8..ddc2a4d 100644 (file)
@@ -115,6 +115,7 @@ c        call chainbuild_cart
       if (nfgtasks.gt.1) then
         call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
       endif
+c      write (iout,*) "itime_mat",itime_mat," imatupdate",imatupdate
       if (mod(itime_mat,imatupdate).eq.0) then
         call make_SCp_inter_list
         call make_SCSC_inter_list
@@ -1476,6 +1477,7 @@ C
      & sigij,r0ij,rcut,sqrij,sss1,sssgrad1
       double precision fcont,fprimcont
       double precision sscale,sscagrad
+      double precision boxshift
 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
 c      do i=iatsc_s,iatsc_e
@@ -1488,6 +1490,7 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
 C Change 12/1/95
         num_conti=0
 C
@@ -1499,9 +1502,13 @@ cd   &                  'iend=',iend(i,iint)
 c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j)) 
             if (itypj.eq.ntyp1) cycle
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
 C Change 12/1/95 to calculate four-body interactions
             rij=xj*xj+yj*yj+zj*zj
             rrij=1.0D0/rij
@@ -1649,6 +1656,7 @@ C
      & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
       logical scheck
       double precision sscale,sscagrad
+      double precision boxshift
 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
 c      do i=iatsc_s,iatsc_e
@@ -1661,6 +1669,7 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
 C
 C Calculate SC interaction energy.
 C
@@ -1668,9 +1677,13 @@ c        do iint=1,nint_gr(i)
 c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             fac_augm=rrij**expon
             e_augm=augm(itypi,itypj)*fac_augm
@@ -1748,6 +1761,7 @@ C
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi,
      & sss1,sssgrad1
       double precision sscale,sscagrad
+      double precision boxshift
 c     double precision rrsave(maxdim)
       logical lprn
       evdw=0.0D0
@@ -1769,6 +1783,7 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -1803,9 +1818,13 @@ c           chip12=0.0D0
 c           alf1=0.0D0
 c           alf2=0.0D0
 c           alf12=0.0D0
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
       include 'COMMON.SPLITELE'
       include 'COMMON.SBRIDGE'
       logical lprn
-      integer xshift,yshift,zshift,subchap
       double precision evdw
       integer itypi,itypj,itypi1,iint,ind,ikont
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
       double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
-     & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
-     & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
+     & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip
       double precision dist,sscale,sscagrad,sscagradlip,sscalelip
+      double precision boxshift
       evdw=0.0D0
 ccccc      energy_dec=.false.
 C      print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -1912,66 +1930,14 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
-C Return atom into box, boxxsize is size of box in x dimension
-c  134   continue
-c        if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
-c        if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
-C Condition for being inside the proper box
-c        if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
-c     &       (xi.lt.((xshift-0.5d0)*boxxsize))) then
-c        go to 134
-c        endif
-c  135   continue
-c        if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
-c        if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
-C Condition for being inside the proper box
-c        if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
-c     &       (yi.lt.((yshift-0.5d0)*boxysize))) then
-c        go to 135
-c        endif
-c  136   continue
-c        if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
-c        if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
-C Condition for being inside the proper box
-c        if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
-c     &       (zi.lt.((zshift-0.5d0)*boxzsize))) then
-c        go to 136
-c        endif
-          xi=mod(xi,boxxsize)
-          if (xi.lt.0) xi=xi+boxxsize
-          yi=mod(yi,boxysize)
-          if (yi.lt.0) yi=yi+boxysize
-          zi=mod(zi,boxzsize)
-          if (zi.lt.0) zi=zi+boxzsize
+        call to_box(xi,yi,zi)
 C define scaling factor for lipids
 
 C        if (positi.le.0) positi=positi+boxzsize
 C        print *,i
 C first for peptide groups
 c for each residue check if it is in lipid or lipid water border area
-       if ((zi.gt.bordlipbot)
-     &.and.(zi.lt.bordliptop)) then
-C the energy transfer exist
-        if (zi.lt.buflipbot) then
-C what fraction I am in
-         fracinbuf=1.0d0-
-     &        ((zi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
-         sslipi=sscalelip(fracinbuf)
-         ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
-        elseif (zi.gt.bufliptop) then
-         fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
-         sslipi=sscalelip(fracinbuf)
-         ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
-        else
-         sslipi=1.0d0
-         ssgradlipi=0.0
-        endif
-       else
-         sslipi=0.0d0
-         ssgradlipi=0.0
-       endif
-
+        call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
 C          xi=xi+xshift*boxxsize
 C          yi=yi+yshift*boxysize
 C          zi=zi+zshift*boxzsize
@@ -1998,15 +1964,15 @@ c              write(iout,*) "PO ZWYKLE", evdwij
               if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') 
      &                        'evdw',i,j,evdwij,' ss'
 C triple bond artifac removal
-             do k=j+1,iend(i,iint) 
+              do k=j+1,iend(i,iint) 
 C search over all next residues
-              if (dyn_ss_mask(k)) then
+                if (dyn_ss_mask(k)) then
 C check if they are cysteins
 C              write(iout,*) 'k=',k
 
 c              write(iout,*) "PRZED TRI", evdwij
-               evdwij_przed_tri=evdwij
-              call triple_ssbond_ene(i,j,k,evdwij)
+                  evdwij_przed_tri=evdwij
+                  call triple_ssbond_ene(i,j,k,evdwij)
 c               if(evdwij_przed_tri.ne.evdwij) then
 c                 write (iout,*) "TRI:", evdwij, evdwij_przed_tri
 c               endif
@@ -2014,30 +1980,30 @@ c               endif
 c              write(iout,*) "PO TRI", evdwij
 C call the energy function that removes the artifical triple disulfide
 C bond the soubroutine is located in ssMD.F
-              evdw=evdw+evdwij             
-              if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+                  evdw=evdw+evdwij             
+                  if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
      &                        'evdw',i,j,evdwij,'tss'
-              endif!dyn_ss_mask(k)
-             enddo! k
+                endif!dyn_ss_mask(k)
+              enddo! k
             ELSE
-            ind=ind+1
-            itypj=iabs(itype(j))
-            if (itypj.eq.ntyp1) cycle
+              ind=ind+1
+              itypj=iabs(itype(j))
+              if (itypj.eq.ntyp1) cycle
 c            dscj_inv=dsc_inv(itypj)
-            dscj_inv=vbld_inv(j+nres)
+              dscj_inv=vbld_inv(j+nres)
 c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
 c     &       1.0d0/vbld(j+nres)
 c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
-            sig0ij=sigma(itypi,itypj)
-            chi1=chi(itypi,itypj)
-            chi2=chi(itypj,itypi)
-            chi12=chi1*chi2
-            chip1=chip(itypi)
-            chip2=chip(itypj)
-            chip12=chip1*chip2
-            alf1=alp(itypi)
-            alf2=alp(itypj)
-            alf12=0.5D0*(alf1+alf2)
+              sig0ij=sigma(itypi,itypj)
+              chi1=chi(itypi,itypj)
+              chi2=chi(itypj,itypi)
+              chi12=chi1*chi2
+              chip1=chip(itypi)
+              chip2=chip(itypj)
+              chip12=chip1*chip2
+              alf1=alp(itypi)
+              alf2=alp(itypj)
+              alf12=0.5D0*(alf1+alf2)
 C For diagnostics only!!!
 c           chi1=0.0D0
 c           chi2=0.0D0
@@ -2048,191 +2014,112 @@ c           chip12=0.0D0
 c           alf1=0.0D0
 c           alf2=0.0D0
 c           alf12=0.0D0
-            xj=c(1,nres+j)
-            yj=c(2,nres+j)
-            zj=c(3,nres+j)
-C Return atom J into box the original box
-c  137   continue
-c        if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
-c        if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
-C Condition for being inside the proper box
-c        if ((xj.gt.((0.5d0)*boxxsize)).or.
-c     &       (xj.lt.((-0.5d0)*boxxsize))) then
-c        go to 137
-c        endif
-c  138   continue
-c        if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
-c        if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
-C Condition for being inside the proper box
-c        if ((yj.gt.((0.5d0)*boxysize)).or.
-c     &       (yj.lt.((-0.5d0)*boxysize))) then
-c        go to 138
-c        endif
-c  139   continue
-c        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
-c        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
-C Condition for being inside the proper box
-c        if ((zj.gt.((0.5d0)*boxzsize)).or.
-c     &       (zj.lt.((-0.5d0)*boxzsize))) then
-c        go to 139
-c        endif
-          xj=mod(xj,boxxsize)
-          if (xj.lt.0) xj=xj+boxxsize
-          yj=mod(yj,boxysize)
-          if (yj.lt.0) yj=yj+boxysize
-          zj=mod(zj,boxzsize)
-          if (zj.lt.0) zj=zj+boxzsize
-       if ((zj.gt.bordlipbot)
-     &.and.(zj.lt.bordliptop)) then
-C the energy transfer exist
-        if (zj.lt.buflipbot) then
-C what fraction I am in
-         fracinbuf=1.0d0-
-     &        ((zj-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
-         sslipj=sscalelip(fracinbuf)
-         ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
-        elseif (zj.gt.bufliptop) then
-         fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
-         sslipj=sscalelip(fracinbuf)
-         ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
-        else
-         sslipj=1.0d0
-         ssgradlipj=0.0
-        endif
-       else
-         sslipj=0.0d0
-         ssgradlipj=0.0
-       endif
-      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
-     &  +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
-      bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
-     &  +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+              xj=c(1,nres+j)
+              yj=c(2,nres+j)
+              zj=c(3,nres+j)
+              call to_box(xj,yj,zj)
+              call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+              aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &          +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+              bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &          +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
 C      write(iout,*) "tu,", i,j,aa_lip(itypi,itypj),bb_lip(itypi,itypj)
 C      if (aa.ne.aa_aq(itypi,itypj)) write(63,'(2e10.5)')
 C     &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj))
 C      if (ssgradlipj.gt.0.0d0) print *,"??WTF??"
 C      print *,sslipi,sslipj,bordlipbot,zi,zj
-      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-      xj_safe=xj
-      yj_safe=yj
-      zj_safe=zj
-      subchap=0
-      do xshift=-1,1
-      do yshift=-1,1
-      do zshift=-1,1
-          xj=xj_safe+xshift*boxxsize
-          yj=yj_safe+yshift*boxysize
-          zj=zj_safe+zshift*boxzsize
-          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-          if(dist_temp.lt.dist_init) then
-            dist_init=dist_temp
-            xj_temp=xj
-            yj_temp=yj
-            zj_temp=zj
-            subchap=1
-          endif
-       enddo
-       enddo
-       enddo
-       if (subchap.eq.1) then
-          xj=xj_temp-xi
-          yj=yj_temp-yi
-          zj=zj_temp-zi
-       else
-          xj=xj_safe-xi
-          yj=yj_safe-yi
-          zj=zj_safe-zi
-       endif
-            dxj=dc_norm(1,nres+j)
-            dyj=dc_norm(2,nres+j)
-            dzj=dc_norm(3,nres+j)
+              xj=boxshift(xj-xi,boxxsize)
+              yj=boxshift(yj-yi,boxysize)
+              zj=boxshift(zj-zi,boxzsize)
+              dxj=dc_norm(1,nres+j)
+              dyj=dc_norm(2,nres+j)
+              dzj=dc_norm(3,nres+j)
 C            xj=xj-xi
 C            yj=yj-yi
 C            zj=zj-zi
 c            write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
 c            write (iout,*) "j",j," dc_norm",
 c     &       dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
-            rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
-            rij=dsqrt(rrij)
-            sss=sscale(1.0d0/rij,r_cut_int)
+              rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+              rij=dsqrt(rrij)
+              sss=sscale(1.0d0/rij,r_cut_int)
 c            write (iout,'(a7,4f8.3)') 
 c    &      "ssscale",sss,((1.0d0/rij)/sigma(itypi,itypj)),r_cut,rlamb
-            if (sss.eq.0.0d0) cycle
-            sssgrad=sscagrad(1.0d0/rij,r_cut_int)
+              if (sss.eq.0.0d0) cycle
+              sssgrad=sscagrad(1.0d0/rij,r_cut_int)
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
-            call sc_angular
-            sigsq=1.0D0/sigsq
-            sig=sig0ij*dsqrt(sigsq)
-            rij_shift=1.0D0/rij-sig+sig0ij
+              call sc_angular
+              sigsq=1.0D0/sigsq
+              sig=sig0ij*dsqrt(sigsq)
+              rij_shift=1.0D0/rij-sig+sig0ij
 c for diagnostics; uncomment
 c            rij_shift=1.2*sig0ij
 C I hate to put IF's in the loops, but here don't have another choice!!!!
-            if (rij_shift.le.0.0D0) then
-              evdw=1.0D20
+              if (rij_shift.le.0.0D0) then
+                evdw=1.0D20
 cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 cd     &        restyp(itypi),i,restyp(itypj),j,
 cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
-              return
-            endif
-            sigder=-sig*sigsq
+                return
+              endif
+              sigder=-sig*sigsq
 c---------------------------------------------------------------
-            rij_shift=1.0D0/rij_shift 
-            fac=rij_shift**expon
+              rij_shift=1.0D0/rij_shift 
+              fac=rij_shift**expon
 C here to start with
 C            if (c(i,3).gt.
-            faclip=fac
-            e1=fac*fac*aa
-            e2=fac*bb
-            evdwij=eps1*eps2rt*eps3rt*(e1+e2)
-            eps2der=evdwij*eps3rt
-            eps3der=evdwij*eps2rt
+              faclip=fac
+              e1=fac*fac*aa
+              e2=fac*bb
+              evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+              eps2der=evdwij*eps3rt
+              eps3der=evdwij*eps2rt
 C       write(63,'(2i3,2e10.3,2f10.5)') i,j,aa,bb, evdwij,
 C     &((sslipi+sslipj)/2.0d0+
 C     &(2.0d0-sslipi-sslipj)/2.0d0)
 c            write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
 c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
-            evdwij=evdwij*eps2rt*eps3rt
-            evdw=evdw+evdwij*sss
-            if (lprn) then
-            sigm=dabs(aa/bb)**(1.0D0/6.0D0)
-            epsi=bb**2/aa
-            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-     &        restyp(itypi),i,restyp(itypj),j,
-     &        epsi,sigm,chi1,chi2,chip1,chip2,
-     &        eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
-     &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
-     &        evdwij
-            endif
+              evdwij=evdwij*eps2rt*eps3rt
+              evdw=evdw+evdwij*sss
+              if (lprn) then
+                sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+                epsi=bb**2/aa
+                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+     &           restyp(itypi),i,restyp(itypj),j,
+     &           epsi,sigm,chi1,chi2,chip1,chip2,
+     &           eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
+     &           om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
+     &           evdwij
+              endif
 
-            if (energy_dec) write (iout,'(a,2i5,3f10.5)') 
+              if (energy_dec) write (iout,'(a,2i5,3f10.5)') 
      &                    'r sss evdw',i,j,1.0d0/rij,sss,evdwij
 
 C Calculate gradient components.
-            e1=e1*eps1*eps2rt**2*eps3rt**2
-            fac=-expon*(e1+evdwij)*rij_shift
-            sigder=fac*sigder
-            fac=rij*fac
+              e1=e1*eps1*eps2rt**2*eps3rt**2
+              fac=-expon*(e1+evdwij)*rij_shift
+              sigder=fac*sigder
+              fac=rij*fac
 c            print '(2i4,6f8.4)',i,j,sss,sssgrad*
 c     &      evdwij,fac,sigma(itypi,itypj),expon
-            fac=fac+evdwij*sssgrad/sss*rij
+              fac=fac+evdwij*sssgrad/sss*rij
 c            fac=0.0d0
 C Calculate the radial part of the gradient
-            gg_lipi(3)=eps1*(eps2rt*eps2rt)
-     &       *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
-     &        (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
-     &       +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
-            gg_lipj(3)=ssgradlipj*gg_lipi(3)
-            gg_lipi(3)=gg_lipi(3)*ssgradlipi
+              gg_lipi(3)=eps1*(eps2rt*eps2rt)
+     &          *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
+     &           (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+     &          +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
+              gg_lipj(3)=ssgradlipj*gg_lipi(3)
+              gg_lipi(3)=gg_lipi(3)*ssgradlipi
 C            gg_lipi(3)=0.0d0
 C            gg_lipj(3)=0.0d0
-            gg(1)=xj*fac
-            gg(2)=yj*fac
-            gg(3)=zj*fac
+              gg(1)=xj*fac
+              gg(2)=yj*fac
+              gg(3)=zj*fac
 C Calculate angular part of the gradient.
 c            call sc_grad_scale(sss)
-            call sc_grad
+              call sc_grad
             ENDIF    ! dyn_ss            
 c          enddo      ! j
 c        enddo        ! iint
@@ -2262,7 +2149,7 @@ C
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
       include 'COMMON.SPLITELE'
-      integer xshift,yshift,zshift,subchap
+      double precision boxshift
       integer icall
       common /srutu/ icall
       logical lprn
@@ -2271,8 +2158,7 @@ C
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
      & xi,yi,zi,fac_augm,e_augm
       double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
-     & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
-     & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip,sssgrad1
+     & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip,sssgrad1
       double precision dist,sscale,sscagrad,sscagradlip,sscalelip
       evdw=0.0D0
 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -2290,41 +2176,14 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
-          xi=mod(xi,boxxsize)
-          if (xi.lt.0) xi=xi+boxxsize
-          yi=mod(yi,boxysize)
-          if (yi.lt.0) yi=yi+boxysize
-          zi=mod(zi,boxzsize)
-          if (zi.lt.0) zi=zi+boxzsize
+        call to_box(xi,yi,zi)
 C define scaling factor for lipids
 
 C        if (positi.le.0) positi=positi+boxzsize
 C        print *,i
 C first for peptide groups
 c for each residue check if it is in lipid or lipid water border area
-       if ((zi.gt.bordlipbot)
-     &.and.(zi.lt.bordliptop)) then
-C the energy transfer exist
-        if (zi.lt.buflipbot) then
-C what fraction I am in
-         fracinbuf=1.0d0-
-     &        ((zi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
-         sslipi=sscalelip(fracinbuf)
-         ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
-        elseif (zi.gt.bufliptop) then
-         fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
-         sslipi=sscalelip(fracinbuf)
-         ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
-        else
-         sslipi=1.0d0
-         ssgradlipi=0.0
-        endif
-       else
-         sslipi=0.0d0
-         ssgradlipi=0.0
-       endif
-
+        call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -2361,131 +2220,77 @@ c           chip12=0.0D0
 c           alf1=0.0D0
 c           alf2=0.0D0
 c           alf12=0.0D0
-C            xj=c(1,nres+j)-xi
-C            yj=c(2,nres+j)-yi
-C            zj=c(3,nres+j)-zi
-          xj=mod(xj,boxxsize)
-          if (xj.lt.0) xj=xj+boxxsize
-          yj=mod(yj,boxysize)
-          if (yj.lt.0) yj=yj+boxysize
-          zj=mod(zj,boxzsize)
-          if (zj.lt.0) zj=zj+boxzsize
-       if ((zj.gt.bordlipbot)
-     &.and.(zj.lt.bordliptop)) then
-C the energy transfer exist
-        if (zj.lt.buflipbot) then
-C what fraction I am in
-         fracinbuf=1.0d0-
-     &        ((zj-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
-         sslipj=sscalelip(fracinbuf)
-         ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
-        elseif (zj.gt.bufliptop) then
-         fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
-         sslipj=sscalelip(fracinbuf)
-         ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
-        else
-         sslipj=1.0d0
-         ssgradlipj=0.0
-        endif
-       else
-         sslipj=0.0d0
-         ssgradlipj=0.0
-       endif
-      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
-     &  +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
-      bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
-     &  +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+           xj=c(1,nres+j)
+           yj=c(2,nres+j)
+           zj=c(3,nres+j)
+           call to_box(xj,yj,zj)
+           call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+           aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+           bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &       +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
 C      if (aa.ne.aa_aq(itypi,itypj)) write(63,'2e10.5') 
 C     &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj))
 C      write(iout,*) "tu,", i,j,aa,bb,aa_lip(itypi,itypj),sslipi,sslipj
-      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-      xj_safe=xj
-      yj_safe=yj
-      zj_safe=zj
-      subchap=0
-      do xshift=-1,1
-      do yshift=-1,1
-      do zshift=-1,1
-          xj=xj_safe+xshift*boxxsize
-          yj=yj_safe+yshift*boxysize
-          zj=zj_safe+zshift*boxzsize
-          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-          if(dist_temp.lt.dist_init) then
-            dist_init=dist_temp
-            xj_temp=xj
-            yj_temp=yj
-            zj_temp=zj
-            subchap=1
-          endif
-       enddo
-       enddo
-       enddo
-       if (subchap.eq.1) then
-          xj=xj_temp-xi
-          yj=yj_temp-yi
-          zj=zj_temp-zi
-       else
-          xj=xj_safe-xi
-          yj=yj_safe-yi
-          zj=zj_safe-zi
-       endif
-            dxj=dc_norm(1,nres+j)
-            dyj=dc_norm(2,nres+j)
-            dzj=dc_norm(3,nres+j)
-            rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
-            rij=dsqrt(rrij)
-            sss=sscale(1.0d0/rij,r_cut_int)
-            if (sss.eq.0.0d0) cycle
-            sssgrad=sscagrad(1.0d0/rij,r_cut_int)
+           xj=boxshift(xj-xi,boxxsize)
+           yj=boxshift(yj-yi,boxysize)
+           zj=boxshift(zj-zi,boxzsize)
+           dxj=dc_norm(1,nres+j)
+           dyj=dc_norm(2,nres+j)
+           dzj=dc_norm(3,nres+j)
+           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+           rij=dsqrt(rrij)
+           sss=sscale(1.0d0/rij,r_cut_int)
+           if (sss.eq.0.0d0) cycle
+           sssgrad=sscagrad(1.0d0/rij,r_cut_int)
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
-            call sc_angular
-            sigsq=1.0D0/sigsq
-            sig=sig0ij*dsqrt(sigsq)
-            rij_shift=1.0D0/rij-sig+r0ij
+           call sc_angular
+           sigsq=1.0D0/sigsq
+           sig=sig0ij*dsqrt(sigsq)
+           rij_shift=1.0D0/rij-sig+r0ij
 C I hate to put IF's in the loops, but here don't have another choice!!!!
-            if (rij_shift.le.0.0D0) then
-              evdw=1.0D20
-              return
-            endif
-            sigder=-sig*sigsq
+           if (rij_shift.le.0.0D0) then
+             evdw=1.0D20
+             return
+           endif
+           sigder=-sig*sigsq
 c---------------------------------------------------------------
-            rij_shift=1.0D0/rij_shift 
-            fac=rij_shift**expon
-            e1=fac*fac*aa
-            e2=fac*bb
-            evdwij=eps1*eps2rt*eps3rt*(e1+e2)
-            eps2der=evdwij*eps3rt
-            eps3der=evdwij*eps2rt
-            fac_augm=rrij**expon
-            e_augm=augm(itypi,itypj)*fac_augm
-            evdwij=evdwij*eps2rt*eps3rt
-            evdw=evdw+evdwij+e_augm
-            if (lprn) then
-            sigm=dabs(aa/bb)**(1.0D0/6.0D0)
-            epsi=bb**2/aa
-            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+           rij_shift=1.0D0/rij_shift 
+           fac=rij_shift**expon
+           e1=fac*fac*aa
+           e2=fac*bb
+           evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+           eps2der=evdwij*eps3rt
+           eps3der=evdwij*eps2rt
+           fac_augm=rrij**expon
+           e_augm=augm(itypi,itypj)*fac_augm
+           evdwij=evdwij*eps2rt*eps3rt
+           evdw=evdw+evdwij+e_augm
+           if (lprn) then
+             sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+             epsi=bb**2/aa
+             write (iout,'(2(a3,i3,2x),17(0pf7.3))')
      &        restyp(itypi),i,restyp(itypj),j,
      &        epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
      &        chi1,chi2,chip1,chip2,
      &        eps1,eps2rt**2,eps3rt**2,
      &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
      &        evdwij+e_augm
-            endif
+           endif
 C Calculate gradient components.
-            e1=e1*eps1*eps2rt**2*eps3rt**2
-            fac=-expon*(e1+evdwij)*rij_shift
-            sigder=fac*sigder
-            fac=rij*fac-2*expon*rrij*e_augm
-            fac=fac+(evdwij+e_augm)*sssgrad/sss*rij
+           e1=e1*eps1*eps2rt**2*eps3rt**2
+           fac=-expon*(e1+evdwij)*rij_shift
+           sigder=fac*sigder
+           fac=rij*fac-2*expon*rrij*e_augm
+           fac=fac+(evdwij+e_augm)*sssgrad/sss*rij
 C Calculate the radial part of the gradient
-            gg(1)=xj*fac
-            gg(2)=yj*fac
-            gg(3)=zj*fac
+           gg(1)=xj*fac
+           gg(2)=yj*fac
+           gg(3)=zj*fac
 C Calculate angular part of the gradient.
 c            call sc_grad_scale(sss)
-            call sc_grad
+           call sc_grad
 c          enddo      ! j
 c        enddo        ! iint
       enddo          ! i
@@ -2636,6 +2441,7 @@ C
       include 'COMMON.IOUNITS'
 c      include 'COMMON.CONTACTS'
       dimension gg(3)
+      double precision boxshift
 cd    print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct
       evdw=0.0D0
 c      do i=iatsc_s,iatsc_e
@@ -2648,6 +2454,7 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
 C
 C Calculate SC interaction energy.
 C
@@ -2657,9 +2464,9 @@ cd   &                  'iend=',iend(i,iint)
 c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=boxshift(c(1,nres+j)-xi,boxxsize)
+            yj=boxshift(c(2,nres+j)-yi,boxysize)
+            zj=boxshift(c(3,nres+j)-zi,boxzsize)
             rij=xj*xj+yj*yj+zj*zj
 c           write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
             r0ij=r0(itypi,itypj)
@@ -2716,7 +2523,7 @@ c      include 'COMMON.CONTACTS'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
       dimension ggg(3)
-      integer xshift,yshift,zshift
+      double precision boxshift
 C      write(iout,*) 'In EELEC_soft_sphere'
       ees=0.0D0
       evdw1=0.0D0
@@ -2732,12 +2539,7 @@ C      write(iout,*) 'In EELEC_soft_sphere'
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
-          xmedi=mod(xmedi,boxxsize)
-          if (xmedi.lt.0) xmedi=xmedi+boxxsize
-          ymedi=mod(ymedi,boxysize)
-          if (ymedi.lt.0) ymedi=ymedi+boxysize
-          zmedi=mod(zmedi,boxzsize)
-          if (zmedi.lt.0) zmedi=zmedi+boxzsize
+        call to_box(xmedi,ymedi,zmedi)
         num_conti=0
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         do j=ielstart(i),ielend(i)
@@ -2754,43 +2556,10 @@ c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
           xj=c(1,j)+0.5D0*dxj
           yj=c(2,j)+0.5D0*dyj
           zj=c(3,j)+0.5D0*dzj
-          xj=mod(xj,boxxsize)
-          if (xj.lt.0) xj=xj+boxxsize
-          yj=mod(yj,boxysize)
-          if (yj.lt.0) yj=yj+boxysize
-          zj=mod(zj,boxzsize)
-          if (zj.lt.0) zj=zj+boxzsize
-      dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
-      xj_safe=xj
-      yj_safe=yj
-      zj_safe=zj
-      isubchap=0
-      do xshift=-1,1
-      do yshift=-1,1
-      do zshift=-1,1
-          xj=xj_safe+xshift*boxxsize
-          yj=yj_safe+yshift*boxysize
-          zj=zj_safe+zshift*boxzsize
-          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-          if(dist_temp.lt.dist_init) then
-            dist_init=dist_temp
-            xj_temp=xj
-            yj_temp=yj
-            zj_temp=zj
-            isubchap=1
-          endif
-       enddo
-       enddo
-       enddo
-       if (isubchap.eq.1) then
-          xj=xj_temp-xmedi
-          yj=yj_temp-ymedi
-          zj=zj_temp-zmedi
-       else
-          xj=xj_safe-xmedi
-          yj=yj_safe-ymedi
-          zj=zj_safe-zmedi
-       endif
+          call to_box(xj,yj,zj)
+          xj=boxshift(xj-xmedi,boxxsize)
+          yj=boxshift(yj-ymedi,boxysize)
+          zj=boxshift(zj-zmedi,boxzsize)
           rij=xj*xj+yj*yj+zj*zj
             sss=sscale(sqrt(rij),r_cut_int)
             sssgrad=sscagrad(sqrt(rij),r_cut_int)
@@ -3786,12 +3555,7 @@ c        end if
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
-          xmedi=mod(xmedi,boxxsize)
-          if (xmedi.lt.0) xmedi=xmedi+boxxsize
-          ymedi=mod(ymedi,boxysize)
-          if (ymedi.lt.0) ymedi=ymedi+boxysize
-          zmedi=mod(zmedi,boxzsize)
-          if (zmedi.lt.0) zmedi=zmedi+boxzsize
+        call to_box(xmedi,ymedi,zmedi)
         num_conti=0
         call eelecij(i,i+2,ees,evdw1,eel_loc)
         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
@@ -3846,13 +3610,7 @@ c        if ((zmedi.gt.((0.5d0)*boxzsize)).or.
 c     &       (zmedi.lt.((-0.5d0)*boxzsize))) then
 c        go to 196
 c        endif
-          xmedi=mod(xmedi,boxxsize)
-          if (xmedi.lt.0) xmedi=xmedi+boxxsize
-          ymedi=mod(ymedi,boxysize)
-          if (ymedi.lt.0) ymedi=ymedi+boxysize
-          zmedi=mod(zmedi,boxzsize)
-          if (zmedi.lt.0) zmedi=zmedi+boxzsize
-
+        call to_box(xmedi,ymedi,zmedi)
 #ifdef FOURBODY
         num_conti=num_cont_hb(i)
 #endif
@@ -3895,12 +3653,7 @@ c     &  .or. itype(i-1).eq.ntyp1
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
-          xmedi=mod(xmedi,boxxsize)
-          if (xmedi.lt.0) xmedi=xmedi+boxxsize
-          ymedi=mod(ymedi,boxysize)
-          if (ymedi.lt.0) ymedi=ymedi+boxysize
-          zmedi=mod(zmedi,boxzsize)
-          if (zmedi.lt.0) zmedi=zmedi+boxzsize
+        call to_box(xmedi,ymedi,zmedi)
 C          xmedi=xmedi+xshift*boxxsize
 C          ymedi=ymedi+yshift*boxysize
 C          zmedi=zmedi+zshift*boxzsize
@@ -3940,7 +3693,7 @@ c        do j=ielstart(i),ielend(i)
 C          do j=16,17
 C          write (iout,*) i,j
 C         if (j.le.1) cycle
-          if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
+        if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
 C changes suggested by Ana to avoid out of bounds
 c     & .or.((j+2).gt.nres)
 c     & .or.((j-1).le.0)
@@ -3948,7 +3701,7 @@ C end of changes by Ana
 c     & .or.itype(j+2).eq.ntyp1
 c     & .or.itype(j-1).eq.ntyp1
      &) cycle
-          call eelecij(i,j,ees,evdw1,eel_loc)
+        call eelecij(i,j,ees,evdw1,eel_loc)
 c        enddo ! j
 #ifdef FOURBODY
         num_cont_hb(i)=num_conti
@@ -4017,10 +3770,9 @@ C-------------------------------------------------------------------------------
      &  ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,
      &  ecosgp,ecosam,ecosbm,ecosgm,ghalf,rlocshield
       double precision a22,a23,a32,a33,geel_loc_ij,geel_loc_ji
-      double precision dist_init,xj_safe,yj_safe,zj_safe,
-     &  xj_temp,yj_temp,zj_temp,dist_temp,xmedi,ymedi,zmedi
+      double precision xmedi,ymedi,zmedi
       double precision sscale,sscagrad,scalar
-
+      double precision boxshift
 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
 #ifdef MOMENT
       double precision scal_el /1.0d0/
@@ -4032,7 +3784,6 @@ C 13-go grudnia roku pamietnego...
       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
      &                   0.0d0,1.0d0,0.0d0,
      &                   0.0d0,0.0d0,1.0d0/
-       integer xshift,yshift,zshift
 c          time00=MPI_Wtime()
 cd      write (iout,*) "eelecij",i,j
 c          ind=ind+1
@@ -4055,73 +3806,12 @@ C          zj=c(3,j)+0.5D0*dzj-zmedi
           xj=c(1,j)+0.5D0*dxj
           yj=c(2,j)+0.5D0*dyj
           zj=c(3,j)+0.5D0*dzj
-          xj=mod(xj,boxxsize)
-          if (xj.lt.0) xj=xj+boxxsize
-          yj=mod(yj,boxysize)
-          if (yj.lt.0) yj=yj+boxysize
-          zj=mod(zj,boxzsize)
-          if (zj.lt.0) zj=zj+boxzsize
-          if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ"
-      dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
-      xj_safe=xj
-      yj_safe=yj
-      zj_safe=zj
-      isubchap=0
-      do xshift=-1,1
-      do yshift=-1,1
-      do zshift=-1,1
-          xj=xj_safe+xshift*boxxsize
-          yj=yj_safe+yshift*boxysize
-          zj=zj_safe+zshift*boxzsize
-          dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
-          if(dist_temp.lt.dist_init) then
-            dist_init=dist_temp
-            xj_temp=xj
-            yj_temp=yj
-            zj_temp=zj
-            isubchap=1
-          endif
-       enddo
-       enddo
-       enddo
-       if (isubchap.eq.1) then
-          xj=xj_temp-xmedi
-          yj=yj_temp-ymedi
-          zj=zj_temp-zmedi
-       else
-          xj=xj_safe-xmedi
-          yj=yj_safe-ymedi
-          zj=zj_safe-zmedi
-       endif
+          call to_box(xj,yj,zj)
+          xj=boxshift(xj-xmedi,boxxsize)
+          yj=boxshift(yj-ymedi,boxysize)
+          zj=boxshift(zj-zmedi,boxzsize)
 C        if ((i+3).lt.j) then !this condition keeps for turn3 and turn4 not subject to PBC
 c  174   continue
-c        if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
-c        if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
-C Condition for being inside the proper box
-c        if ((xj.gt.((0.5d0)*boxxsize)).or.
-c     &       (xj.lt.((-0.5d0)*boxxsize))) then
-c        go to 174
-c        endif
-c  175   continue
-c        if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
-c        if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
-C Condition for being inside the proper box
-c        if ((yj.gt.((0.5d0)*boxysize)).or.
-c     &       (yj.lt.((-0.5d0)*boxysize))) then
-c        go to 175
-c        endif
-c  176   continue
-c        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
-c        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
-C Condition for being inside the proper box
-c        if ((zj.gt.((0.5d0)*boxzsize)).or.
-c     &       (zj.lt.((-0.5d0)*boxzsize))) then
-c        go to 176
-c        endif
-C        endif !endPBC condintion
-C        xj=xj-xmedi
-C        yj=yj-ymedi
-C        zj=zj-zmedi
           rij=xj*xj+yj*yj+zj*zj
 
           sss=sscale(dsqrt(rij),r_cut_int)
@@ -5569,7 +5259,7 @@ C
       include 'COMMON.IOUNITS'
       include 'COMMON.CONTROL'
       dimension ggg(3)
-      integer xshift,yshift,zshift
+      double precision boxshift
       evdw2=0.0D0
       evdw2_14=0.0d0
       r0_scp=4.5d0
@@ -5612,12 +5302,7 @@ c        if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
 c     &       (zi.lt.((zshift-0.5d0)*boxzsize))) then
 c        go to 136
 c        endif
-          xi=mod(xi,boxxsize)
-          if (xi.lt.0) xi=xi+boxxsize
-          yi=mod(yi,boxysize)
-          if (yi.lt.0) yi=yi+boxysize
-          zi=mod(zi,boxzsize)
-          if (zi.lt.0) zi=zi+boxzsize
+          call to_box(xi,yi,zi)
 C          xi=xi+xshift*boxxsize
 C          yi=yi+yshift*boxysize
 C          zi=zi+zshift*boxzsize
@@ -5634,67 +5319,10 @@ C Uncomment following three lines for Ca-p interactions
           xj=c(1,j)
           yj=c(2,j)
           zj=c(3,j)
-c  174   continue
-c        if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
-c        if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
-C Condition for being inside the proper box
-c        if ((xj.gt.((0.5d0)*boxxsize)).or.
-c     &       (xj.lt.((-0.5d0)*boxxsize))) then
-c        go to 174
-c        endif
-c  175   continue
-c        if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
-c        if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
-cC Condition for being inside the proper box
-c        if ((yj.gt.((0.5d0)*boxysize)).or.
-c     &       (yj.lt.((-0.5d0)*boxysize))) then
-c        go to 175
-c        endif
-c  176   continue
-c        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
-c        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
-C Condition for being inside the proper box
-c        if ((zj.gt.((0.5d0)*boxzsize)).or.
-c     &       (zj.lt.((-0.5d0)*boxzsize))) then
-c        go to 176
-          xj=mod(xj,boxxsize)
-          if (xj.lt.0) xj=xj+boxxsize
-          yj=mod(yj,boxysize)
-          if (yj.lt.0) yj=yj+boxysize
-          zj=mod(zj,boxzsize)
-          if (zj.lt.0) zj=zj+boxzsize
-      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-      xj_safe=xj
-      yj_safe=yj
-      zj_safe=zj
-      subchap=0
-      do xshift=-1,1
-      do yshift=-1,1
-      do zshift=-1,1
-          xj=xj_safe+xshift*boxxsize
-          yj=yj_safe+yshift*boxysize
-          zj=zj_safe+zshift*boxzsize
-          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-          if(dist_temp.lt.dist_init) then
-            dist_init=dist_temp
-            xj_temp=xj
-            yj_temp=yj
-            zj_temp=zj
-            subchap=1
-          endif
-       enddo
-       enddo
-       enddo
-       if (subchap.eq.1) then
-          xj=xj_temp-xi
-          yj=yj_temp-yi
-          zj=zj_temp-zi
-       else
-          xj=xj_safe-xi
-          yj=yj_safe-yi
-          zj=zj_safe-zi
-       endif
-c c       endif
+          call to_box(xj,yj,zj)
+          xj=boxshift(xj-xi,boxxsize)
+          yj=boxshift(yj-yi,boxysize)
+          zj=boxshift(zj-zi,boxzsize)
 C          xj=xj-xi
 C          yj=yj-yi
 C          zj=zj-zi
@@ -5716,32 +5344,6 @@ C
           ggg(1)=xj*fac
           ggg(2)=yj*fac
           ggg(3)=zj*fac
-cgrad          if (j.lt.i) then
-cd          write (iout,*) 'j<i'
-C Uncomment following three lines for SC-p interactions
-c           do k=1,3
-c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
-c           enddo
-cgrad          else
-cd          write (iout,*) 'j>i'
-cgrad            do k=1,3
-cgrad              ggg(k)=-ggg(k)
-C Uncomment following line for SC-p interactions
-c             gradx_scp(k,j)=gradx_scp(k,j)-ggg(k)
-cgrad            enddo
-cgrad          endif
-cgrad          do k=1,3
-cgrad            gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k)
-cgrad          enddo
-cgrad          kstart=min0(i+1,j)
-cgrad          kend=max0(i-1,j-1)
-cd        write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend
-cd        write (iout,*) ggg(1),ggg(2),ggg(3)
-cgrad          do k=kstart,kend
-cgrad            do l=1,3
-cgrad              gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l)
-cgrad            enddo
-cgrad          enddo
           do k=1,3
             gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
             gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
       include 'COMMON.IOUNITS'
       include 'COMMON.CONTROL'
       include 'COMMON.SPLITELE'
-      integer xshift,yshift,zshift
       double precision ggg(3)
       integer i,iint,j,k,iteli,itypj,subchap,ikont
       double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
      & fac,e1,e2,rij
       double precision evdw2,evdw2_14,evdwij
-      double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
-     & dist_temp, dist_init
       double precision sscale,sscagrad
+      double precision boxshift
       evdw2=0.0D0
       evdw2_14=0.0d0
 c        print *,boxxsize,boxysize,boxzsize,'wymiary pudla'
@@ -5801,43 +5401,7 @@ c      do i=iatscp_s,iatscp_e
         xi=0.5D0*(c(1,i)+c(1,i+1))
         yi=0.5D0*(c(2,i)+c(2,i+1))
         zi=0.5D0*(c(3,i)+c(3,i+1))
-          xi=mod(xi,boxxsize)
-          if (xi.lt.0) xi=xi+boxxsize
-          yi=mod(yi,boxysize)
-          if (yi.lt.0) yi=yi+boxysize
-          zi=mod(zi,boxzsize)
-          if (zi.lt.0) zi=zi+boxzsize
-c          xi=xi+xshift*boxxsize
-c          yi=yi+yshift*boxysize
-c          zi=zi+zshift*boxzsize
-c        print *,xi,yi,zi,'polozenie i'
-C Return atom into box, boxxsize is size of box in x dimension
-c  134   continue
-c        if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
-c        if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
-C Condition for being inside the proper box
-c        if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
-c     &       (xi.lt.((xshift-0.5d0)*boxxsize))) then
-c        go to 134
-c        endif
-c  135   continue
-c          print *,xi,boxxsize,"pierwszy"
-
-c        if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
-c        if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
-C Condition for being inside the proper box
-c        if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
-c     &       (yi.lt.((yshift-0.5d0)*boxysize))) then
-c        go to 135
-c        endif
-c  136   continue
-c        if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
-c        if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
-C Condition for being inside the proper box
-c        if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
-c     &       (zi.lt.((zshift-0.5d0)*boxzsize))) then
-c        go to 136
-c        endif
+        call to_box(xi,yi,zi)
 c        do iint=1,nscp_gr(i)
 
 c        do j=iscpstart(i,iint),iscpend(i,iint)
@@ -5851,68 +5415,10 @@ C Uncomment following three lines for Ca-p interactions
           xj=c(1,j)
           yj=c(2,j)
           zj=c(3,j)
-          xj=mod(xj,boxxsize)
-          if (xj.lt.0) xj=xj+boxxsize
-          yj=mod(yj,boxysize)
-          if (yj.lt.0) yj=yj+boxysize
-          zj=mod(zj,boxzsize)
-          if (zj.lt.0) zj=zj+boxzsize
-c  174   continue
-c        if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
-c        if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
-C Condition for being inside the proper box
-c        if ((xj.gt.((0.5d0)*boxxsize)).or.
-c     &       (xj.lt.((-0.5d0)*boxxsize))) then
-c        go to 174
-c        endif
-c  175   continue
-c        if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
-c        if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
-cC Condition for being inside the proper box
-c        if ((yj.gt.((0.5d0)*boxysize)).or.
-c     &       (yj.lt.((-0.5d0)*boxysize))) then
-c        go to 175
-c        endif
-c  176   continue
-c        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
-c        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
-C Condition for being inside the proper box
-c        if ((zj.gt.((0.5d0)*boxzsize)).or.
-c     &       (zj.lt.((-0.5d0)*boxzsize))) then
-c        go to 176
-c        endif
-CHERE IS THE CALCULATION WHICH MIRROR IMAGE IS THE CLOSEST ONE
-      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-      xj_safe=xj
-      yj_safe=yj
-      zj_safe=zj
-      subchap=0
-      do xshift=-1,1
-      do yshift=-1,1
-      do zshift=-1,1
-          xj=xj_safe+xshift*boxxsize
-          yj=yj_safe+yshift*boxysize
-          zj=zj_safe+zshift*boxzsize
-          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-          if(dist_temp.lt.dist_init) then
-            dist_init=dist_temp
-            xj_temp=xj
-            yj_temp=yj
-            zj_temp=zj
-            subchap=1
-          endif
-       enddo
-       enddo
-       enddo
-       if (subchap.eq.1) then
-          xj=xj_temp-xi
-          yj=yj_temp-yi
-          zj=zj_temp-zi
-       else
-          xj=xj_safe-xi
-          yj=yj_safe-yi
-          zj=zj_safe-zi
-       endif
+          call to_box(xj,yj,zj)
+          xj=boxshift(xj-xi,boxxsize)
+          yj=boxshift(yj-yi,boxysize)
+          zj=boxshift(zj-zi,boxzsize)
 c          print *,xj,yj,zj,'polozenie j'
           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 c          print *,rrij
@@ -13638,3 +13144,103 @@ C-----------------------------------------------------------------------
       endif
       return
       end
+c------------------------------------------------------------------------
+      double precision function boxshift(x,boxsize)
+      implicit none
+      double precision x,boxsize
+      double precision xtemp
+      xtemp=dmod(x,boxsize)
+      if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then
+        boxshift=xtemp-boxsize
+      else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then
+        boxshift=xtemp+boxsize
+      else
+        boxshift=xtemp
+      endif
+      return
+      end
+c--------------------------------------------------------------------------
+      subroutine closest_img(xi,yi,zi,xj,yj,zj)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      integer xshift,yshift,zshift,subchap
+      double precision dist_init,xj_safe,yj_safe,zj_safe,
+     & xj_temp,yj_temp,zj_temp,dist_temp
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+      subchap=0
+      do xshift=-1,1
+        do yshift=-1,1
+          do zshift=-1,1
+            xj=xj_safe+xshift*boxxsize
+            yj=yj_safe+yshift*boxysize
+            zj=zj_safe+zshift*boxzsize
+            dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+            if(dist_temp.lt.dist_init) then
+              dist_init=dist_temp
+              xj_temp=xj
+              yj_temp=yj
+              zj_temp=zj
+              subchap=1
+            endif
+          enddo
+        enddo
+      enddo
+      if (subchap.eq.1) then
+        xj=xj_temp-xi
+        yj=yj_temp-yi
+        zj=zj_temp-zi
+      else
+        xj=xj_safe-xi
+        yj=yj_safe-yi
+        zj=zj_safe-zi
+      endif
+      return
+      end
+c--------------------------------------------------------------------------
+      subroutine to_box(xi,yi,zi)
+      implicit none
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      double precision xi,yi,zi
+      xi=dmod(xi,boxxsize)
+      if (xi.lt.0.0d0) xi=xi+boxxsize
+      yi=dmod(yi,boxysize)
+      if (yi.lt.0.0d0) yi=yi+boxysize
+      zi=dmod(zi,boxzsize)
+      if (zi.lt.0.0d0) zi=zi+boxzsize
+      return
+      end
+c--------------------------------------------------------------------------
+      subroutine lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+      implicit none
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      double precision xi,yi,zi,sslipi,ssgradlipi
+      double precision fracinbuf
+      double precision sscalelip,sscagradlip
+
+      if ((zi.gt.bordlipbot).and.(zi.lt.bordliptop)) then
+C the energy transfer exist
+        if (zi.lt.buflipbot) then
+C what fraction I am in
+          fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+          sslipi=sscalelip(fracinbuf)
+          ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+        elseif (zi.gt.bufliptop) then
+          fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
+          sslipi=sscalelip(fracinbuf)
+          ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+        else
+          sslipi=1.0d0
+          ssgradlipi=0.0
+        endif
+      else
+        sslipi=0.0d0
+        ssgradlipi=0.0
+      endif
+      return
+      end
index 3e662cc..557f435 100644 (file)
@@ -845,11 +845,12 @@ c     overlapping residues left, or false otherwise (success)
       integer ioverlap(maxres),ioverlap_last
       data redfac /0.5D0/
 
-      write (iout,*) "overlap_sc_list"
+c      write (iout,*) "overlap_sc_list"
       ioverlap_last=0
 C Check for SC-SC overlaps and mark residues
 c      print *,'>>overlap_sc nnt=',nnt,' nct=',nct
       ind=0
+c      write(iout,*) "iatsc_s",iatsc_s," iatsc_e",iatsc_e
       do i=iatsc_s,iatsc_e
         itypi=iabs(itype(i))
         itypi1=iabs(itype(i+1))
@@ -884,7 +885,7 @@ c            write (iout,*) "i,j",i,j," itypi,itypj",itypi,itypj
           else 
            rcomp=sigma(itypi,itypj)
           endif
-c         print '(2(a3,2i3),a3,2f10.5)',
+c          write (iout,'(2(a3,2i3),a3,2f10.5)'),
 c     &        ' i=',i,iti,' j=',j,itj,' d=',dist(nres+i,nres+j)
 c     &        ,rcomp
             xj=c(1,nres+j)-xi
@@ -896,17 +897,23 @@ c     &        ,rcomp
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
             call sc_angular
+c            write (iout,*) "dxj",dxj," dyj",dyj," dzj",dzj
+c            write (iout,*) "erij",erij
+c            write (iout,*) "om1",om1," om2",om2," om12",om12,
+c     &       " faceps1",faceps1," eps1",eps1
+c            write (iout,*) "sigsq",sigsq
             sigsq=1.0D0/sigsq
             sig=sig0ij*dsqrt(sigsq)
             rij_shift=1.0D0/rij-sig+sig0ij
-
-ct          if ( 1.0/rij .lt. redfac*rcomp .or. 
-ct     &       rij_shift.le.0.0D0 ) then
-c           write (iout,'(a,i3,a,i3,a,f10.5,a,3f10.5)')
+c            write (iout,*) "rij_shift",rij_shift
+c          if ( 1.0/rij .lt. redfac*rcomp .or. 
+c     &       rij_shift.le.0.0D0 ) then
+c          write (iout,'(a,i3,a,i3,a,f10.5,a,3f10.5)')
 c     &     'overlap SC-SC: i=',i,' j=',j,
 c     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
 c     &     rcomp,1.0/rij,rij_shift
             if ( rij_shift.le.0.0D0 ) then
+c                    write (iout,*) "overlap",i,j
           ioverlap_last=ioverlap_last+1
           ioverlap(ioverlap_last)=i         
           do k=1,ioverlap_last-1
diff --git a/source/unres/src-HCD-5D/int_from_cart.F b/source/unres/src-HCD-5D/int_from_cart.F
new file mode 100644 (file)
index 0000000..2ae703b
--- /dev/null
@@ -0,0 +1,269 @@
+      subroutine int_from_cart(lside,lprn)
+      implicit none
+      include 'DIMENSIONS'
+#ifdef MPI
+      include "mpif.h"
+#endif 
+      include 'COMMON.LOCAL'
+      include 'COMMON.VAR'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.GEO'
+      include 'COMMON.NAMES'
+      include 'COMMON.CONTROL'
+      include 'COMMON.SETUP'
+      double precision dist,alpha,beta
+      character*3 seq,atom,res
+      character*80 card
+      double precision sccor(3,50)
+      integer rescode
+      logical lside,lprn
+      integer i,j,iti
+      double precision di,cosfac2,sinfac2,cosfac,sinfac
+#ifdef MPI
+      if(me.eq.king.or..not.out1file)then
+#endif
+       if (lprn) then 
+        write (iout,'(/a)') 
+     &  'Internal coordinates calculated from crystal structure.'
+        if (lside) then 
+          write (iout,'(8a)') '  Res  ','       dvb','     Theta',
+     & '       Phi','    Dsc_id','       Dsc','     Alpha',
+     & '     Omega'
+        else 
+          write (iout,'(4a)') '  Res  ','       dvb','     Theta',
+     & '       Phi'
+        endif
+       endif
+#ifdef MPI
+      endif
+#endif
+      do i=1,nres-1
+        iti=itype(i)
+        if (iti.ne.ntyp1 .and. itype(i+1).ne.ntyp1 .and. 
+     &      (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0)) then
+          write (iout,'(a,i4)') 'Bad Cartesians for residue',i
+ctest          stop
+        endif
+        vbld(i+1)=dist(i,i+1)
+        vbld_inv(i+1)=1.0d0/vbld(i+1)
+c        write (iout,*) "i",i+1," vbld",vbld(i+1)," vbld_inv",
+c     &      vbld_inv(i+1)
+        if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
+        if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
+      enddo
+c      if (unres_pdb) then
+c        if (itype(1).eq.21) then
+c          theta(3)=90.0d0*deg2rad
+c          phi(4)=180.0d0*deg2rad
+c          vbld(2)=3.8d0
+c          vbld_inv(2)=1.0d0/vbld(2)
+c        endif
+c        if (itype(nres).eq.21) then
+c          theta(nres)=90.0d0*deg2rad
+c          phi(nres)=180.0d0*deg2rad
+c          vbld(nres)=3.8d0
+c          vbld_inv(nres)=1.0d0/vbld(2)
+c        endif
+c      endif
+c      print *,"A TU2"
+      if (lside) then
+        do i=2,nres-1
+          do j=1,3
+            c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
+     &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
+          enddo
+          iti=itype(i)
+          di=dist(i,nres+i)
+          vbld(i+nres)=di
+          if (itype(i).ne.10) then
+            vbld_inv(i+nres)=1.0d0/di
+          else
+            vbld_inv(i+nres)=0.0d0
+          endif
+          if (iti.ne.10) then
+            alph(i)=alpha(nres+i,i,maxres2)
+            omeg(i)=beta(nres+i,i,maxres2,i+1)
+          endif
+          if(me.eq.king.or..not.out1file)then
+           if (lprn)
+     &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
+     &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
+     &     rad2deg*alph(i),rad2deg*omeg(i)
+          endif
+        enddo
+        if (lprn) then
+           i=nres
+           iti=itype(nres)
+           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
+     &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
+     &     rad2deg*alph(i),rad2deg*omeg(i)
+        endif
+      else if (lprn) then
+        do i=2,nres
+          iti=itype(i)
+          if(me.eq.king.or..not.out1file)
+     &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
+     &     rad2deg*theta(i),rad2deg*phi(i)
+        enddo
+      endif
+      return
+      end
+c-------------------------------------------------------------------------------
+      subroutine sc_loc_geom(lprn)
+      implicit none
+      include 'DIMENSIONS'
+#ifdef MPI
+      include "mpif.h"
+#endif 
+      include 'COMMON.LOCAL'
+      include 'COMMON.VAR'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.GEO'
+      include 'COMMON.NAMES'
+      include 'COMMON.CONTROL'
+      include 'COMMON.SETUP'
+      double precision x_prime(3),y_prime(3),z_prime(3)
+      logical lprn
+      integer i,j,it
+      double precision xx,yy,zz,cosfac,cosfac2,sinfac,sinfac2
+      do i=1,nres-1
+        do j=1,3
+          dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
+        enddo
+c        write (iout,*) "i",i," dc",(dc_norm(j,i),j=1,3),
+c     &    " vbld",vbld_inv(i+1)
+      enddo
+      do i=2,nres-1
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+          do j=1,3
+            dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
+          enddo
+c          write (iout,*) "i",i," dc",(dc_norm(j,i+nres),j=1,3),
+c     &    " vbld",vbld_inv(i+nres)
+        else
+          do j=1,3
+            dc_norm(j,i+nres)=0.0d0
+          enddo
+        endif
+      enddo
+      do i=2,nres-1
+        costtab(i+1) =dcos(theta(i+1))
+        sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
+        cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
+        sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
+        cosfac2=0.5d0/(1.0d0+costtab(i+1))
+        cosfac=dsqrt(cosfac2)
+        sinfac2=0.5d0/(1.0d0-costtab(i+1))
+        sinfac=dsqrt(sinfac2)
+        it=itype(i)
+c        write (iout,*) "i",i," costab",costtab(i+1),
+c     &    " sintab",sinttab(i+1)
+c        write (iout,*) "dc_norm_b",(dc_norm(j,i-1),j=1,3)
+c        write (iout,*) "dc_norm_s",(dc_norm(j,i+nres),j=1,3)
+        if (it.ne.10 .and. itype(i).ne.ntyp1) then
+c
+C  Compute the axes of tghe local cartesian coordinates system; store in
+c   x_prime, y_prime and z_prime 
+c
+        do j=1,3
+          x_prime(j) = 0.00
+          y_prime(j) = 0.00
+          z_prime(j) = 0.00
+        enddo
+        do j = 1,3
+          x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
+          y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
+        enddo
+c        write (iout,*) "x_prime",(x_prime(j),j=1,3)
+c        write (iout,*) "y_prime",(y_prime(j),j=1,3)
+        call vecpr(x_prime,y_prime,z_prime)
+c        write (iout,*) "z_prime",(z_prime(j),j=1,3)
+c
+C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
+C to local coordinate system. Store in xx, yy, zz.
+c
+        xx=0.0d0
+        yy=0.0d0
+        zz=0.0d0
+        do j = 1,3
+          xx = xx + x_prime(j)*dc_norm(j,i+nres)
+          yy = yy + y_prime(j)*dc_norm(j,i+nres)
+          zz = zz + z_prime(j)*dc_norm(j,i+nres)
+        enddo
+
+        xxref(i)=xx
+        yyref(i)=yy
+        zzref(i)=zz
+        else
+        xxref(i)=0.0d0
+        yyref(i)=0.0d0
+        zzref(i)=0.0d0
+        endif
+      enddo
+      if (lprn) then
+#ifdef MPI
+        if (me.eq.king.or..not.out1file) then
+#endif
+        write (iout,*) "xxref,yyref,zzref"
+        do i=2,nres
+          write (iout,'(a3,i4,3f10.5)') 
+     &     restyp(itype(i)),i,xxref(i),yyref(i),zzref(i)
+        enddo
+#ifdef MPI
+        endif
+#endif
+      endif
+      return
+      end
+c---------------------------------------------------------------------------
+      subroutine sccenter(ires,nscat,sccor)
+      implicit none
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      integer i,j,ires,nscat
+      double precision sccor(3,50)
+      double precision sccmj
+      do j=1,3
+        sccmj=0.0D0
+        do i=1,nscat
+          sccmj=sccmj+sccor(j,i) 
+        enddo
+        dc(j,ires)=sccmj/nscat
+      enddo
+      return
+      end
+c---------------------------------------------------------------------------
+      subroutine bond_regular
+      implicit none
+      include 'DIMENSIONS'   
+      include 'COMMON.VAR'
+      include 'COMMON.LOCAL'      
+      include 'COMMON.INTERACT'
+      include 'COMMON.CHAIN'
+      integer i,i1,i2
+      do i=1,nres-1
+       vbld(i+1)=vbl
+       vbld_inv(i+1)=vblinv
+       vbld(i+1+nres)=dsc(iabs(itype(i+1)))
+       vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1)))
+c       print *,vbld(i+1),vbld(i+1+nres)
+      enddo
+c Adam 2/26/20 Alter virtual bonds for non-blocking end groups of each chain
+      do i=1,nchain
+        i1=chain_border(1,i)
+        i2=chain_border(2,i)
+        if (i1.gt.1) then
+          vbld(i1)=vbld(i1)/2
+          vbld_inv(i1)=vbld_inv(i1)*2
+        endif
+        if (i2.lt.nres) then
+          vbld(i2+1)=vbld(i2+1)/2
+          vbld_inv(i2+1)=vbld_inv(i2+1)*2
+        endif
+      enddo
+      return
+      end
index a56e4f8..cea54c4 100644 (file)
@@ -635,8 +635,7 @@ c     v(25)=4.0D0
         endif
       enddo
       nvarx=k
-      write (iout,*) "Variables set up nvarx",nvarx
-      write (iout,*) "Before energy minimization"
+      call chainbuild_cart
       call etotal(energia(0))
       call enerprint(energia(0))
 #ifdef LBFGS
diff --git a/source/unres/src-HCD-5D/readpdb-mult.F b/source/unres/src-HCD-5D/readpdb-mult.F
new file mode 100644 (file)
index 0000000..8fd1da8
--- /dev/null
@@ -0,0 +1,635 @@
+      subroutine readpdb
+C Read the PDB file and convert the peptide geometry into virtual-chain 
+C geometry.
+      implicit none
+      include 'DIMENSIONS'
+      include 'COMMON.LOCAL'
+      include 'COMMON.VAR'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.GEO'
+      include 'COMMON.NAMES'
+      include 'COMMON.CONTROL'
+      include 'COMMON.FRAG'
+      include 'COMMON.SETUP'
+      include 'COMMON.SBRIDGE'
+      character*3 seq,atom,res
+      character*80 card
+      double precision sccor(3,50)
+      double precision e1(3),e2(3),e3(3)
+      integer rescode,iterter(maxres),cou
+      logical fail
+      integer i,j,iii,ires,ires_old,ishift,ishift1,ibeg
+      double precision dcj
+      bfac=0.0d0
+      do i=1,maxres
+         iterter(i)=0
+      enddo
+      ibeg=1
+      ishift1=0
+      lsecondary=.false.
+      nhfrag=0
+      nbfrag=0
+      do
+        read (ipdbin,'(a80)',end=10) card
+        if (card(:5).eq.'HELIX') then
+         nhfrag=nhfrag+1
+         lsecondary=.true.
+         read(card(22:25),*) hfrag(1,nhfrag)
+         read(card(34:37),*) hfrag(2,nhfrag)
+        endif
+        if (card(:5).eq.'SHEET') then
+         nbfrag=nbfrag+1
+         lsecondary=.true.
+         read(card(24:26),*) bfrag(1,nbfrag)
+         read(card(35:37),*) bfrag(2,nbfrag)
+crc----------------------------------------
+crc  to be corrected !!!
+         bfrag(3,nbfrag)=bfrag(1,nbfrag)
+         bfrag(4,nbfrag)=bfrag(2,nbfrag)
+crc----------------------------------------
+        endif
+        if (card(:3).eq.'END') then
+          goto 10
+        else if (card(:3).eq.'TER') then
+C End current chain
+          ires_old=ires+2
+          itype(ires_old-1)=ntyp1 
+          iterter(ires_old-1)=1
+          itype(ires_old)=ntyp1
+          iterter(ires_old)=1
+          ibeg=2
+          write (iout,*) "Chain ended",ires,ishift,ires_old
+          if (unres_pdb) then
+            do j=1,3
+              dc(j,ires)=sccor(j,iii)
+            enddo
+          else 
+            call sccenter(ires,iii,sccor)
+          endif
+        endif
+C Fish out the ATOM cards.
+c        if (index(card(1:4),'ATOM').gt.0) then  
+c          read (card(14:16),'(a3)') atom
+c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
+C Calculate the CM of the preceding residue.
+        read (card(23:26),*) ires
+        read (card(18:20),'(a3)') res
+        if (ires-ishift+ishift1.ne.ires_old) then
+          if (ibeg.eq.0) then
+            if (unres_pdb) then
+              do j=1,3
+                dc(j,ires+nres)=sccor(j,iii)
+              enddo
+            else
+              call sccenter(ires,iii,sccor)
+            endif
+          endif
+C Start new residue.
+c            write (iout,'(a80)') card
+          if (ibeg.eq.1) then
+            ishift=ires-1
+            if (res.ne.'GLY' .and. res.ne. 'ACE') then
+              ishift=ishift-1
+              itype(1)=ntyp1
+            endif
+c             write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
+            ibeg=0          
+          else if (ibeg.eq.2) then
+c Start a new chain
+            ishift=-ires_old+ires-1
+c              write (iout,*) "New chain started",ires,ishift
+            ibeg=0
+          endif
+        else
+          ires=ires-ishift
+c            write (2,*) "ires",ires," ishift",ish
+        endif
+        if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
+     &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
+          if (res.eq.'ACE') then
+            itype(ires)=10
+          else
+            itype(ires)=rescode(ires,res,0)
+          endif
+          read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+          read(card(61:66),*) bfac(ires)
+c            if(me.eq.king.or..not.out1file)
+c     &       write (iout,'(2i3,2x,a,3f8.3)') 
+c     &       ires,itype(ires),res,(c(j,ires),j=1,3)
+          iii=1
+          do j=1,3
+            sccor(j,iii)=c(j,ires)
+          enddo
+        else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and. 
+     &             atom.ne.'N  ' .and. atom.ne.'C   ') then
+            iii=iii+1
+          read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+        endif
+      enddo
+   10 if(me.eq.king.or..not.out1file) 
+     & write (iout,'(a,i5)') ' Nres: ',ires
+C Calculate dummy residue coordinates inside the "chain" of a multichain
+C system
+      nres=ires
+      do i=2,nres-1
+c        write (iout,*) i,itype(i),itype(i+1),ntyp1,iterter(i)
+        if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
+         if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
+C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
+C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
+C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
+           if (unres_pdb) then
+C 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
+            print *,i,'a tu?'
+            do j=1,3
+             c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
+            enddo
+           else   !unres_pdb
+           do j=1,3
+             dcj=(c(j,i-2)-c(j,i-3))/2.0
+            if (dcj.eq.0) dcj=1.23591524223
+             c(j,i)=c(j,i-1)+dcj
+             c(j,nres+i)=c(j,i)
+           enddo     
+          endif   !unres_pdb
+         else     !itype(i+1).eq.ntyp1
+          if (unres_pdb) then
+C 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
+            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+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 !itype(i+1).eq.ntyp1
+        endif  !itype.eq.ntyp1
+      enddo
+      write (iout,*) "After loop in readpbd"
+C Calculate the CM of the last side chain.
+      if (unres_pdb) then
+        do j=1,3
+          dc(j,ires)=sccor(j,iii)
+        enddo
+      else 
+        call sccenter(ires,iii,sccor)
+      endif
+      nsup=nres
+      nstart_sup=1
+      if (itype(nres).ne.10) then
+        nres=nres+1
+        itype(nres)=ntyp1
+        if (unres_pdb) then
+C 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*(-e1(j)+e2(j))/sqrt(2.0d0)
+          enddo
+        else
+        do j=1,3
+          dcj=(c(j,nres-2)-c(j,nres-3))/2.0
+            if (dcj.eq.0) dcj=1.23591524223
+          c(j,nres)=c(j,nres-1)+dcj
+          c(j,2*nres)=c(j,nres)
+        enddo
+        endif
+      endif
+      do i=2,nres-1
+        do j=1,3
+          c(j,i+nres)=dc(j,i)
+        enddo
+      enddo
+      do j=1,3
+        c(j,nres+1)=c(j,1)
+        c(j,2*nres)=c(j,nres)
+      enddo
+      if (itype(1).eq.ntyp1) then
+        nsup=nsup-1
+        nstart_sup=2
+        if (unres_pdb) then
+C 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*(e1(j)-e2(j))/dsqrt(2.0d0)
+          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
+C Calculate internal coordinates.
+      if(me.eq.king.or..not.out1file)then
+       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),
+     &    (c(j,nres+ires),j=1,3)
+       enddo
+      endif
+      call flush(iout)
+c      write(iout,*)"before int_from_cart nres",nres
+      call int_from_cart(.true.,.false.)
+      do i=1,nres
+        thetaref(i)=theta(i)
+        phiref(i)=phi(i)
+      enddo
+      do i=1,nres-1
+        do j=1,3
+          dc(j,i)=c(j,i+1)-c(j,i)
+          dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+        enddo
+c        write (iout,*) i,(dc(j,i),j=1,3),(dc_norm(j,i),j=1,3),
+c     &   vbld_inv(i+1)
+      enddo
+      do i=2,nres-1
+        do j=1,3
+          dc(j,i+nres)=c(j,i+nres)-c(j,i)
+          dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+        enddo
+c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
+c     &   vbld_inv(i+nres)
+      enddo
+      call sc_loc_geom(.false.)
+      call int_from_cart1(.false.)
+c      call chainbuild
+C Copy the coordinates to reference coordinates
+      do i=1,nres
+        do j=1,3
+          cref(j,i)=c(j,i)
+          cref(j,i+nres)=c(j,i+nres)
+        enddo
+      enddo
+  100 format (//'              alpha-carbon coordinates       ',
+     &          '     centroid coordinates'/
+     1          '       ', 6X,'X',11X,'Y',11X,'Z',
+     &                          10X,'X',11X,'Y',11X,'Z')
+  110 format (a,'(',i3,')',6f12.5)
+cc enddiag
+      do j=1,nbfrag     
+        do i=1,4                                                       
+         bfrag(i,j)=bfrag(i,j)-ishift
+        enddo
+      enddo
+
+      do j=1,nhfrag
+        do i=1,2
+         hfrag(i,j)=hfrag(i,j)-ishift
+        enddo
+      enddo
+      return
+      end
+c---------------------------------------------------------------------------
+      subroutine readpdb_template(k)
+C Read the PDB file for read_constr_homology with read2sigma
+C and convert the peptide geometry into virtual-chain geometry.
+      implicit none
+      include 'DIMENSIONS'
+      include 'COMMON.LOCAL'
+      include 'COMMON.VAR'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.GEO'
+      include 'COMMON.NAMES'
+      include 'COMMON.CONTROL'
+      include 'COMMON.FRAG'
+      include 'COMMON.SETUP'
+      integer i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
+     &  ishift_pdb,ires_ca
+      logical lprn /.false./,fail
+      double precision e1(3),e2(3),e3(3)
+      double precision dcj,efree_temp
+      character*3 seq,res
+      character*5 atom
+      character*80 card
+      double precision sccor(3,20)
+      integer rescode,iterter(maxres)
+      do i=1,maxres
+         iterter(i)=0
+      enddo
+      ibeg=1
+      ishift1=0
+      ishift=0
+c      write (2,*) "UNRES_PDB",unres_pdb
+      ires=0
+      ires_old=0
+      iii=0
+      lsecondary=.false.
+      nhfrag=0
+      nbfrag=0
+      do
+        read (ipdbin,'(a80)',end=10) card
+        if (card(:3).eq.'END') then
+          goto 10
+        else if (card(:3).eq.'TER') then
+C End current chain
+          ires_old=ires+2
+          itype(ires_old-1)=ntyp1 
+          iterter(ires_old-1)=1
+          itype(ires_old)=ntyp1
+          iterter(ires_old)=1
+          ibeg=2
+c          write (iout,*) "Chain ended",ires,ishift,ires_old
+          if (unres_pdb) then
+            do j=1,3
+              dc(j,ires)=sccor(j,iii)
+            enddo
+          else 
+            call sccenter(ires,iii,sccor)
+          endif
+        endif
+C Fish out the ATOM cards.
+        if (index(card(1:4),'ATOM').gt.0) then  
+          read (card(12:16),*) atom
+c          write (iout,*) "! ",atom," !",ires
+c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
+          read (card(23:26),*) ires
+          read (card(18:20),'(a3)') res
+c          write (iout,*) "ires",ires,ires-ishift+ishift1,
+c     &      " ires_old",ires_old
+c          write (iout,*) "ishift",ishift," ishift1",ishift1
+c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
+          if (ires-ishift+ishift1.ne.ires_old) then
+C Calculate the CM of the preceding residue.
+            if (ibeg.eq.0) then
+              if (unres_pdb) then
+                do j=1,3
+                  dc(j,ires)=sccor(j,iii)
+                enddo
+              else
+                call sccenter(ires_old,iii,sccor)
+              endif
+              iii=0
+            endif
+C Start new residue.
+            if (res.eq.'Cl-' .or. res.eq.'Na+') then
+              ires=ires_old
+              cycle
+            else if (ibeg.eq.1) then
+c              write (iout,*) "BEG ires",ires
+              ishift=ires-1
+              if (res.ne.'GLY' .and. res.ne. 'ACE') then
+                ishift=ishift-1
+                itype(1)=ntyp1
+              endif
+              ires=ires-ishift+ishift1
+              ires_old=ires
+c              write (iout,*) "ishift",ishift," ires",ires,
+c     &         " ires_old",ires_old
+c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
+              ibeg=0          
+            else if (ibeg.eq.2) then
+c Start a new chain
+              ishift=-ires_old+ires-1
+              ires=ires_old+1
+c              write (iout,*) "New chain started",ires,ishift
+              ibeg=0          
+            else
+              ishift=ishift-(ires-ishift+ishift1-ires_old-1)
+              ires=ires-ishift+ishift1
+              ires_old=ires
+            endif
+            if (res.eq.'ACE' .or. res.eq.'NHE') then
+              itype(ires)=10
+            else
+              itype(ires)=rescode(ires,res,0)
+            endif
+          else
+            ires=ires-ishift+ishift1
+          endif
+c          write (iout,*) "ires_old",ires_old," ires",ires
+c          if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+c            ishift1=ishift1+1
+c          endif
+c          write (2,*) "ires",ires," res ",res," ity",ity
+          if (atom.eq.'CA' .or. atom.eq.'CH3' .or. 
+     &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
+            read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+c            write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
+#ifdef DEBUG
+            write (iout,'(2i3,2x,a,3f8.3)') 
+     &      ires,itype(ires),res,(c(j,ires),j=1,3)
+#endif
+            iii=iii+1
+            do j=1,3
+              sccor(j,iii)=c(j,ires)
+            enddo
+            if (ishift.ne.0) then
+              ires_ca=ires+ishift-ishift1
+            else
+              ires_ca=ires
+            endif
+c            write (*,*) card(23:27),ires,itype(ires)
+          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.
+     &             atom.ne.'OXT' .and. atom(:2).ne.'3H') then
+c            write (iout,*) "sidechain ",atom
+            iii=iii+1
+            read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+          endif
+        endif
+      enddo
+   10 if(me.eq.king.or..not.out1file) 
+     & write (iout,'(a,i5)') ' Nres: ',ires
+C Calculate dummy residue coordinates inside the "chain" of a multichain
+C system
+      nres=ires
+      do i=2,nres-1
+c        write (iout,*) i,itype(i),itype(i+1)
+        if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
+         if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
+C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
+C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
+C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
+           if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the last dummy residue
+            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
+            do j=1,3
+             c(j,i)=c(j,i-1)-1.9d0*e2(j)
+            enddo
+           else   !unres_pdb
+           do j=1,3
+             dcj=(c(j,i-2)-c(j,i-3))/2.0
+            if (dcj.eq.0) dcj=1.23591524223
+             c(j,i)=c(j,i-1)+dcj
+             c(j,nres+i)=c(j,i)
+           enddo     
+          endif   !unres_pdb
+         else     !itype(i+1).eq.ntyp1
+          if (unres_pdb) then
+C 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
+            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+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 !itype(i+1).eq.ntyp1
+        endif  !itype.eq.ntyp1
+      enddo
+C Calculate the CM of the last side chain.
+      if (unres_pdb) then
+        do j=1,3
+          dc(j,ires)=sccor(j,iii)
+        enddo
+      else
+        call sccenter(ires,iii,sccor)
+      endif
+      nsup=nres
+      nstart_sup=1
+      if (itype(nres).ne.10) then
+        nres=nres+1
+        itype(nres)=ntyp1
+        if (unres_pdb) then
+C 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
+        do j=1,3
+          dcj=(c(j,nres-2)-c(j,nres-3))/2.0
+            if (dcj.eq.0) dcj=1.23591524223
+          c(j,nres)=c(j,nres-1)+dcj
+          c(j,2*nres)=c(j,nres)
+        enddo
+      endif
+      endif
+      do i=2,nres-1
+        do j=1,3
+          c(j,i+nres)=dc(j,i)
+        enddo
+      enddo
+      do j=1,3
+        c(j,nres+1)=c(j,1)
+        c(j,2*nres)=c(j,nres)
+      enddo
+      if (itype(1).eq.ntyp1) then
+        nsup=nsup-1
+        nstart_sup=2
+        if (unres_pdb) then
+C 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
+C Copy the coordinates to reference coordinates
+c      do i=1,2*nres
+c        do j=1,3
+c          cref(j,i)=c(j,i)
+c        enddo
+c      enddo
+C Calculate internal coordinates.
+      if (out_template_coord) then
+      write (iout,'(/a)') 
+     &  "Cartesian coordinates of the reference structure"
+      write (iout,'(a,3(3x,a5),5x,3(3x,a5))') 
+     & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+      do ires=1,nres
+        write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') 
+     &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
+     &    (c(j,ires+nres),j=1,3)
+      enddo
+      endif
+C Calculate internal coordinates.
+      call int_from_cart(.true.,.true.)
+      call sc_loc_geom(.false.)
+      do i=1,nres
+        thetaref(i)=theta(i)
+        phiref(i)=phi(i)
+      enddo
+      do i=1,nres-1
+        do j=1,3
+          dc(j,i)=c(j,i+1)-c(j,i)
+          dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+        enddo
+      enddo
+      do i=2,nres-1
+        do j=1,3
+          dc(j,i+nres)=c(j,i+nres)-c(j,i)
+          dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+        enddo
+c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
+c     &   vbld_inv(i+nres)
+      enddo
+      do i=1,nres
+        do j=1,3
+          cref(j,i)=c(j,i)
+          cref(j,i+nres)=c(j,i+nres)
+        enddo
+      enddo
+      do i=1,2*nres
+        do j=1,3
+          chomo(j,i,k)=c(j,i)
+        enddo
+      enddo
+
+      return
+      end
+      
index 943d67d..9c8462d 100644 (file)
@@ -1,7 +1,7 @@
       subroutine readpdb
 C Read the PDB file and convert the peptide geometry into virtual-chain 
 C geometry.
-      implicit none
+      implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'COMMON.LOCAL'
       include 'COMMON.VAR'
@@ -11,27 +11,33 @@ C geometry.
       include 'COMMON.GEO'
       include 'COMMON.NAMES'
       include 'COMMON.CONTROL'
-      include 'COMMON.FRAG'
+      include 'COMMON.DISTFIT'
       include 'COMMON.SETUP'
-      include 'COMMON.SBRIDGE'
-      character*3 seq,atom,res
+      include 'COMMON.FRAG'
+      integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
+     &  ishift_pdb
+      logical lprn /.false./,fail
+      double precision e1(3),e2(3),e3(3)
+      double precision dcj,efree_temp
+      character*3 seq,res
+      character*5 atom
       character*80 card
       double precision sccor(3,50)
-      double precision e1(3),e2(3),e3(3)
-      integer rescode,iterter(maxres),cou
-      logical fail
-      integer i,j,iii,ires,ires_old,ishift,ibeg
-      double precision dcj
-      bfac=0.0d0
-      do i=1,maxres
-         iterter(i)=0
-      enddo
+      integer rescode
+      efree_temp=0.0d0
       ibeg=1
+      ishift1=0
+      ishift=0
+c      write (2,*) "UNRES_PDB",unres_pdb
+      ires=0
+      ires_old=0
+      iii=0
       lsecondary=.false.
       nhfrag=0
       nbfrag=0
       do
         read (ipdbin,'(a80)',end=10) card
+c        write (iout,'(a)') card
         if (card(:5).eq.'HELIX') then
          nhfrag=nhfrag+1
          lsecondary=.true.
@@ -49,145 +55,112 @@ crc  to be corrected !!!
          bfrag(4,nbfrag)=bfrag(2,nbfrag)
 crc----------------------------------------
         endif
-        if (card(:3).eq.'END') then
-          goto 10
-        else if (card(:3).eq.'TER') then
-C End current chain
-          ires_old=ires+2
-          itype(ires_old-1)=ntyp1 
-          iterter(ires_old-1)=1
-          itype(ires_old)=ntyp1
-          iterter(ires_old)=1
-          ibeg=2
-          write (iout,*) "Chain ended",ires,ishift,ires_old
-          if (unres_pdb) then
-            do j=1,3
-              dc(j,ires)=sccor(j,iii)
-            enddo
-          else 
-            call sccenter(ires,iii,sccor)
-          endif
-        endif
+        if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
+c Read free energy
+        if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
 C Fish out the ATOM cards.
         if (index(card(1:4),'ATOM').gt.0) then  
-          read (card(14:16),'(a3)') atom
-          if (atom.eq.'CA' .or. atom.eq.'CH3') then
+          read (card(12:16),*) atom
+c          write (iout,*) "! ",atom," !",ires
+c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
+          read (card(23:26),*) ires
+          read (card(18:20),'(a3)') res
+c          write (iout,*) "ires",ires,ires-ishift+ishift1,
+c     &      " ires_old",ires_old
+c          write (iout,*) "ishift",ishift," ishift1",ishift1
+c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
+          if (ires-ishift+ishift1.ne.ires_old) then
 C Calculate the CM of the preceding residue.
+c            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
             if (ibeg.eq.0) then
+c              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)=sccor(j,iii)
                 enddo
               else
-                call sccenter(ires,iii,sccor)
+                call sccenter(ires_old,iii,sccor)
               endif
+              iii=0
             endif
 C Start new residue.
-c            write (iout,'(a80)') card
-            read (card(23:26),*) ires
-            read (card(18:20),'(a3)') res
-            if (ibeg.eq.1) then
+            if (res.eq.'Cl-' .or. res.eq.'Na+') then
+              ires=ires_old
+              cycle
+            else if (ibeg.eq.1) then
+c              write (iout,*) "BEG ires",ires
               ishift=ires-1
               if (res.ne.'GLY' .and. res.ne. 'ACE') then
                 ishift=ishift-1
                 itype(1)=ntyp1
               endif
-c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
+              ires=ires-ishift+ishift1
+              ires_old=ires
+c              write (iout,*) "ishift",ishift," ires",ires,
+c     &         " ires_old",ires_old
               ibeg=0          
-            else if (ibeg.eq.2) then
-c Start a new chain
-              ishift=-ires_old+ires-1
-c              write (iout,*) "New chain started",ires,ishift
-              ibeg=0
+            else
+              ishift=ishift-(ires-ishift+ishift1-ires_old-1)
+              ires=ires-ishift+ishift1
+              ires_old=ires
             endif
-            ires=ires-ishift
-c            write (2,*) "ires",ires," ishift",ishift
-            if (res.eq.'ACE') then
+            if (res.eq.'ACE' .or. res.eq.'NHE') then
               itype(ires)=10
             else
               itype(ires)=rescode(ires,res,0)
             endif
+          else
+            ires=ires-ishift+ishift1
+          endif
+c          write (iout,*) "ires_old",ires_old," ires",ires
+          if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+c            ishift1=ishift1+1
+          endif
+c          write (2,*) "ires",ires," res ",res," ity",ity
+          if (atom.eq.'CA' .or. atom.eq.'CH3' .or. 
+     &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
-            read(card(61:66),*) bfac(ires)
-c            if(me.eq.king.or..not.out1file)
-c     &       write (iout,'(2i3,2x,a,3f8.3)') 
-c     &       ires,itype(ires),res,(c(j,ires),j=1,3)
-            iii=1
+c            write (iout,*) "backbone ",atom 
+#ifdef DEBUG
+            write (iout,'(2i3,2x,a,3f8.3)') 
+     &      ires,itype(ires),res,(c(j,ires),j=1,3)
+#endif
+            iii=iii+1
             do j=1,3
               sccor(j,iii)=c(j,ires)
             enddo
-          else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and. 
-     &             atom.ne.'N  ' .and. atom.ne.'C   ') then
+            if (ishift.ne.0) then
+              ires_ca=ires+ishift-ishift1
+            else
+              ires_ca=ires
+            endif
+c            write (*,*) card(23:27),ires,itype(ires)
+          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.
+     &             atom.ne.'OXT' .and. atom(:2).ne.'3H') then
+c            write (iout,*) "sidechain ",atom
             iii=iii+1
             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
           endif
         endif
       enddo
-   10 if(me.eq.king.or..not.out1file) 
-     & write (iout,'(a,i5)') ' Nres: ',ires
-C Calculate dummy residue coordinates inside the "chain" of a multichain
-C system
-      nres=ires
-      do i=2,nres-1
-c        write (iout,*) i,itype(i),itype(i+1),ntyp1,iterter(i)
-        if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
-         if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
-C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
-C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
-C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
-           if (unres_pdb) then
-C 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
-            print *,i,'a tu?'
-            do j=1,3
-             c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
-            enddo
-           else   !unres_pdb
-           do j=1,3
-             dcj=(c(j,i-2)-c(j,i-3))/2.0
-            if (dcj.eq.0) dcj=1.23591524223
-             c(j,i)=c(j,i-1)+dcj
-             c(j,nres+i)=c(j,i)
-           enddo     
-          endif   !unres_pdb
-         else     !itype(i+1).eq.ntyp1
-          if (unres_pdb) then
-C 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
-            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+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 !itype(i+1).eq.ntyp1
-        endif  !itype.eq.ntyp1
-      enddo
-      write (iout,*) "After loop in readpbd"
+   10 continue
+#ifdef DEBUG
+      write (iout,'(a,i5)') ' Number of residues found: ',ires
+#endif
+      if (ires.eq.0) return
 C Calculate the CM of the last side chain.
+      if (iii.gt.0)  then
       if (unres_pdb) then
         do j=1,3
           dc(j,ires)=sccor(j,iii)
         enddo
-      else 
+      else
         call sccenter(ires,iii,sccor)
       endif
+      endif
+      nres=ires
       nsup=nres
       nstart_sup=1
       if (itype(nres).ne.10) then
@@ -202,12 +175,11 @@ C 2/15/2013 by Adam: corrected insertion of the last dummy residue
             e2(3)=0.0d0
           endif
           do j=1,3
-            c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
+            c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
           enddo
         else
         do j=1,3
-          dcj=(c(j,nres-2)-c(j,nres-3))/2.0
-            if (dcj.eq.0) dcj=1.23591524223
+          dcj=c(j,nres-2)-c(j,nres-3)
           c(j,nres)=c(j,nres-1)+dcj
           c(j,2*nres)=c(j,nres)
         enddo
@@ -234,27 +206,46 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue
             e2(3)=0.0d0
           endif
           do j=1,3
-            c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0)
+            c(j,1)=c(j,2)-3.8d0*e2(j)
           enddo
         else
         do j=1,3
-          dcj=(c(j,4)-c(j,3))/2.0
+          dcj=c(j,4)-c(j,3)
           c(j,1)=c(j,2)-dcj
           c(j,nres+1)=c(j,1)
         enddo
         endif
       endif
+C Copy the coordinates to reference coordinates
+c      do i=1,2*nres
+c        do j=1,3
+c          cref(j,i)=c(j,i)
+c        enddo
+c      enddo
+C Calculate internal coordinates.
+      if (lprn) then
+      write (iout,'(/a)') 
+     &  "Cartesian coordinates of the reference structure"
+      write (iout,'(a,3(3x,a5),5x,3(3x,a5))') 
+     & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+      do ires=1,nres
+        write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') 
+     &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
+     &    (c(j,ires+nres),j=1,3)
+      enddo
+      endif
 C Calculate internal coordinates.
       if(me.eq.king.or..not.out1file)then
+       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)') 
      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
      &    (c(j,nres+ires),j=1,3)
        enddo
       endif
-      call flush(iout)
-c      write(iout,*)"before int_from_cart nres",nres
       call int_from_cart(.true.,.false.)
+      call sc_loc_geom(.false.)
       do i=1,nres
         thetaref(i)=theta(i)
         phiref(i)=phi(i)
@@ -264,8 +255,6 @@ c      write(iout,*)"before int_from_cart nres",nres
           dc(j,i)=c(j,i+1)-c(j,i)
           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
         enddo
-c        write (iout,*) i,(dc(j,i),j=1,3),(dc_norm(j,i),j=1,3),
-c     &   vbld_inv(i+1)
       enddo
       do i=2,nres-1
         do j=1,3
@@ -275,22 +264,15 @@ c     &   vbld_inv(i+1)
 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 c     &   vbld_inv(i+nres)
       enddo
-      call sc_loc_geom(.false.)
-      call int_from_cart1(.false.)
 c      call chainbuild
 C Copy the coordinates to reference coordinates
-      do i=1,nres
+      do i=1,2*nres
         do j=1,3
           cref(j,i)=c(j,i)
-          cref(j,i+nres)=c(j,i+nres)
         enddo
       enddo
-  100 format (//'              alpha-carbon coordinates       ',
-     &          '     centroid coordinates'/
-     1          '       ', 6X,'X',11X,'Y',11X,'Z',
-     &                          10X,'X',11X,'Y',11X,'Z')
-  110 format (a,'(',i3,')',6f12.5)
-cc enddiag
+
+
       do j=1,nbfrag     
         do i=1,4                                                       
          bfrag(i,j)=bfrag(i,j)-ishift
@@ -302,283 +284,14 @@ cc enddiag
          hfrag(i,j)=hfrag(i,j)-ishift
         enddo
       enddo
+      ishift_pdb=ishift
       return
       end
-c---------------------------------------------------------------------------
-      subroutine int_from_cart(lside,lprn)
-      implicit none
-      include 'DIMENSIONS'
-#ifdef MPI
-      include "mpif.h"
-#endif 
-      include 'COMMON.LOCAL'
-      include 'COMMON.VAR'
-      include 'COMMON.CHAIN'
-      include 'COMMON.INTERACT'
-      include 'COMMON.IOUNITS'
-      include 'COMMON.GEO'
-      include 'COMMON.NAMES'
-      include 'COMMON.CONTROL'
-      include 'COMMON.SETUP'
-      double precision dist,alpha,beta
-      character*3 seq,atom,res
-      character*80 card
-      double precision sccor(3,50)
-      integer rescode
-      logical lside,lprn
-      integer i,j,iti
-      double precision di,cosfac2,sinfac2,cosfac,sinfac
-#ifdef MPI
-      if(me.eq.king.or..not.out1file)then
-#endif
-       if (lprn) then 
-        write (iout,'(/a)') 
-     &  'Internal coordinates calculated from crystal structure.'
-        if (lside) then 
-          write (iout,'(8a)') '  Res  ','       dvb','     Theta',
-     & '       Phi','    Dsc_id','       Dsc','     Alpha',
-     & '     Omega'
-        else 
-          write (iout,'(4a)') '  Res  ','       dvb','     Theta',
-     & '       Phi'
-        endif
-       endif
-#ifdef MPI
-      endif
-#endif
-      do i=1,nres-1
-        iti=itype(i)
-        if (iti.ne.ntyp1 .and. itype(i+1).ne.ntyp1 .and. 
-     &      (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0)) then
-          write (iout,'(a,i4)') 'Bad Cartesians for residue',i
-ctest          stop
-        endif
-        vbld(i+1)=dist(i,i+1)
-        vbld_inv(i+1)=1.0d0/vbld(i+1)
-c        write (iout,*) "i",i+1," vbld",vbld(i+1)," vbld_inv",
-c     &      vbld_inv(i+1)
-        if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
-        if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
-      enddo
-c      if (unres_pdb) then
-c        if (itype(1).eq.21) then
-c          theta(3)=90.0d0*deg2rad
-c          phi(4)=180.0d0*deg2rad
-c          vbld(2)=3.8d0
-c          vbld_inv(2)=1.0d0/vbld(2)
-c        endif
-c        if (itype(nres).eq.21) then
-c          theta(nres)=90.0d0*deg2rad
-c          phi(nres)=180.0d0*deg2rad
-c          vbld(nres)=3.8d0
-c          vbld_inv(nres)=1.0d0/vbld(2)
-c        endif
-c      endif
-c      print *,"A TU2"
-      if (lside) then
-        do i=2,nres-1
-          do j=1,3
-            c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
-     &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
-          enddo
-          iti=itype(i)
-          di=dist(i,nres+i)
-          vbld(i+nres)=di
-          if (itype(i).ne.10) then
-            vbld_inv(i+nres)=1.0d0/di
-          else
-            vbld_inv(i+nres)=0.0d0
-          endif
-          if (iti.ne.10) then
-            alph(i)=alpha(nres+i,i,maxres2)
-            omeg(i)=beta(nres+i,i,maxres2,i+1)
-          endif
-          if(me.eq.king.or..not.out1file)then
-           if (lprn)
-     &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
-     &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
-     &     rad2deg*alph(i),rad2deg*omeg(i)
-          endif
-        enddo
-        if (lprn) then
-           i=nres
-           iti=itype(nres)
-           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
-     &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
-     &     rad2deg*alph(i),rad2deg*omeg(i)
-        endif
-      else if (lprn) then
-        do i=2,nres
-          iti=itype(i)
-          if(me.eq.king.or..not.out1file)
-     &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
-     &     rad2deg*theta(i),rad2deg*phi(i)
-        enddo
-      endif
-      return
-      end
-c-------------------------------------------------------------------------------
-      subroutine sc_loc_geom(lprn)
-      implicit none
-      include 'DIMENSIONS'
-#ifdef MPI
-      include "mpif.h"
-#endif 
-      include 'COMMON.LOCAL'
-      include 'COMMON.VAR'
-      include 'COMMON.CHAIN'
-      include 'COMMON.INTERACT'
-      include 'COMMON.IOUNITS'
-      include 'COMMON.GEO'
-      include 'COMMON.NAMES'
-      include 'COMMON.CONTROL'
-      include 'COMMON.SETUP'
-      double precision x_prime(3),y_prime(3),z_prime(3)
-      logical lprn
-      integer i,j,it
-      double precision xx,yy,zz,cosfac,cosfac2,sinfac,sinfac2
-      do i=1,nres-1
-        do j=1,3
-          dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
-        enddo
-c        write (iout,*) "i",i," dc",(dc_norm(j,i),j=1,3),
-c     &    " vbld",vbld_inv(i+1)
-      enddo
-      do i=2,nres-1
-        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
-          do j=1,3
-            dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
-          enddo
-c          write (iout,*) "i",i," dc",(dc_norm(j,i+nres),j=1,3),
-c     &    " vbld",vbld_inv(i+nres)
-        else
-          do j=1,3
-            dc_norm(j,i+nres)=0.0d0
-          enddo
-        endif
-      enddo
-      do i=2,nres-1
-        costtab(i+1) =dcos(theta(i+1))
-        sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
-        cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
-        sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
-        cosfac2=0.5d0/(1.0d0+costtab(i+1))
-        cosfac=dsqrt(cosfac2)
-        sinfac2=0.5d0/(1.0d0-costtab(i+1))
-        sinfac=dsqrt(sinfac2)
-        it=itype(i)
-c        write (iout,*) "i",i," costab",costtab(i+1),
-c     &    " sintab",sinttab(i+1)
-c        write (iout,*) "dc_norm_b",(dc_norm(j,i-1),j=1,3)
-c        write (iout,*) "dc_norm_s",(dc_norm(j,i+nres),j=1,3)
-        if (it.ne.10 .and. itype(i).ne.ntyp1) then
-c
-C  Compute the axes of tghe local cartesian coordinates system; store in
-c   x_prime, y_prime and z_prime 
-c
-        do j=1,3
-          x_prime(j) = 0.00
-          y_prime(j) = 0.00
-          z_prime(j) = 0.00
-        enddo
-        do j = 1,3
-          x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
-          y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
-        enddo
-c        write (iout,*) "x_prime",(x_prime(j),j=1,3)
-c        write (iout,*) "y_prime",(y_prime(j),j=1,3)
-        call vecpr(x_prime,y_prime,z_prime)
-c        write (iout,*) "z_prime",(z_prime(j),j=1,3)
-c
-C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
-C to local coordinate system. Store in xx, yy, zz.
-c
-        xx=0.0d0
-        yy=0.0d0
-        zz=0.0d0
-        do j = 1,3
-          xx = xx + x_prime(j)*dc_norm(j,i+nres)
-          yy = yy + y_prime(j)*dc_norm(j,i+nres)
-          zz = zz + z_prime(j)*dc_norm(j,i+nres)
-        enddo
-
-        xxref(i)=xx
-        yyref(i)=yy
-        zzref(i)=zz
-        else
-        xxref(i)=0.0d0
-        yyref(i)=0.0d0
-        zzref(i)=0.0d0
-        endif
-      enddo
-      if (lprn) then
-#ifdef MPI
-        if (me.eq.king.or..not.out1file) then
-#endif
-        write (iout,*) "xxref,yyref,zzref"
-        do i=2,nres
-          write (iout,'(a3,i4,3f10.5)') 
-     &     restyp(itype(i)),i,xxref(i),yyref(i),zzref(i)
-        enddo
-#ifdef MPI
-        endif
-#endif
-      endif
-      return
-      end
-c---------------------------------------------------------------------------
-      subroutine sccenter(ires,nscat,sccor)
-      implicit none
-      include 'DIMENSIONS'
-      include 'COMMON.CHAIN'
-      integer i,j,ires,nscat
-      double precision sccor(3,50)
-      double precision sccmj
-      do j=1,3
-        sccmj=0.0D0
-        do i=1,nscat
-          sccmj=sccmj+sccor(j,i) 
-        enddo
-        dc(j,ires)=sccmj/nscat
-      enddo
-      return
-      end
-c---------------------------------------------------------------------------
-      subroutine bond_regular
-      implicit none
-      include 'DIMENSIONS'   
-      include 'COMMON.VAR'
-      include 'COMMON.LOCAL'      
-      include 'COMMON.INTERACT'
-      include 'COMMON.CHAIN'
-      integer i,i1,i2
-      do i=1,nres-1
-       vbld(i+1)=vbl
-       vbld_inv(i+1)=vblinv
-       vbld(i+1+nres)=dsc(iabs(itype(i+1)))
-       vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1)))
-c       print *,vbld(i+1),vbld(i+1+nres)
-      enddo
-c Adam 2/26/20 Alter virtual bonds for non-blocking end groups of each chain
-      do i=1,nchain
-        i1=chain_border(1,i)
-        i2=chain_border(2,i)
-        if (i1.gt.1) then
-          vbld(i1)=vbld(i1)/2
-          vbld_inv(i1)=vbld_inv(i1)*2
-        endif
-        if (i2.lt.nres) then
-          vbld(i2+1)=vbld(i2+1)/2
-          vbld_inv(i2+1)=vbld_inv(i2+1)*2
-        endif
-      enddo
-      return
-      end
-c---------------------------------------------------------------------------
+c-----------------------------------------------------------------------
       subroutine readpdb_template(k)
-C Read the PDB file for read_constr_homology with read2sigma
+C Read the PDB file with gaps for read_constr_homology with read2sigma
 C and convert the peptide geometry into virtual-chain geometry.
-      implicit none
+      implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'COMMON.LOCAL'
       include 'COMMON.VAR'
@@ -588,21 +301,20 @@ C and convert the peptide geometry into virtual-chain geometry.
       include 'COMMON.GEO'
       include 'COMMON.NAMES'
       include 'COMMON.CONTROL'
-      include 'COMMON.FRAG'
+      include 'COMMON.DISTFIT'
       include 'COMMON.SETUP'
-      integer i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
-     &  ishift_pdb,ires_ca
+      include 'COMMON.MD'
+      integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
+     &  ishift_pdb
       logical lprn /.false./,fail
       double precision e1(3),e2(3),e3(3)
       double precision dcj,efree_temp
       character*3 seq,res
       character*5 atom
       character*80 card
-      double precision sccor(3,20)
-      integer rescode,iterter(maxres)
-      do i=1,maxres
-         iterter(i)=0
-      enddo
+      double precision sccor(3,50)
+      integer rescode
+      efree_temp=0.0d0
       ibeg=1
       ishift1=0
       ishift=0
@@ -613,27 +325,10 @@ c      write (2,*) "UNRES_PDB",unres_pdb
       lsecondary=.false.
       nhfrag=0
       nbfrag=0
-      do
+      do 
         read (ipdbin,'(a80)',end=10) card
-        if (card(:3).eq.'END') then
-          goto 10
-        else if (card(:3).eq.'TER') then
-C End current chain
-          ires_old=ires+2
-          itype(ires_old-1)=ntyp1 
-          iterter(ires_old-1)=1
-          itype(ires_old)=ntyp1
-          iterter(ires_old)=1
-          ibeg=2
-c          write (iout,*) "Chain ended",ires,ishift,ires_old
-          if (unres_pdb) then
-            do j=1,3
-              dc(j,ires)=sccor(j,iii)
-            enddo
-          else 
-            call sccenter(ires,iii,sccor)
-          endif
-        endif
+c        write (iout,'(a)') card
+        if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
 C Fish out the ATOM cards.
         if (index(card(1:4),'ATOM').gt.0) then  
           read (card(12:16),*) atom
@@ -647,7 +342,9 @@ c          write (iout,*) "ishift",ishift," ishift1",ishift1
 c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
           if (ires-ishift+ishift1.ne.ires_old) then
 C Calculate the CM of the preceding residue.
+c            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
             if (ibeg.eq.0) then
+c              write (iout,*) "Calculating sidechain center iii",iii
               if (unres_pdb) then
                 do j=1,3
                   dc(j,ires)=sccor(j,iii)
@@ -672,13 +369,6 @@ c              write (iout,*) "BEG ires",ires
               ires_old=ires
 c              write (iout,*) "ishift",ishift," ires",ires,
 c     &         " ires_old",ires_old
-c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
-              ibeg=0          
-            else if (ibeg.eq.2) then
-c Start a new chain
-              ishift=-ires_old+ires-1
-              ires=ires_old+1
-c              write (iout,*) "New chain started",ires,ishift
               ibeg=0          
             else
               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
@@ -694,14 +384,14 @@ c              write (iout,*) "New chain started",ires,ishift
             ires=ires-ishift+ishift1
           endif
 c          write (iout,*) "ires_old",ires_old," ires",ires
-c          if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+          if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 c            ishift1=ishift1+1
-c          endif
+          endif
 c          write (2,*) "ires",ires," res ",res," ity",ity
           if (atom.eq.'CA' .or. atom.eq.'CH3' .or. 
      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
-c            write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
+c            write (iout,*) "backbone ",atom 
 #ifdef DEBUG
             write (iout,'(2i3,2x,a,3f8.3)') 
      &      ires,itype(ires),res,(c(j,ires),j=1,3)
@@ -726,61 +416,13 @@ c            write (iout,*) "sidechain ",atom
           endif
         endif
       enddo
-   10 if(me.eq.king.or..not.out1file) 
-     & write (iout,'(a,i5)') ' Nres: ',ires
-C Calculate dummy residue coordinates inside the "chain" of a multichain
-C system
-      nres=ires
-      do i=2,nres-1
-c        write (iout,*) i,itype(i),itype(i+1)
-        if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
-         if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
-C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
-C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
-C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
-           if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the last dummy residue
-            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
-            do j=1,3
-             c(j,i)=c(j,i-1)-1.9d0*e2(j)
-            enddo
-           else   !unres_pdb
-           do j=1,3
-             dcj=(c(j,i-2)-c(j,i-3))/2.0
-            if (dcj.eq.0) dcj=1.23591524223
-             c(j,i)=c(j,i-1)+dcj
-             c(j,nres+i)=c(j,i)
-           enddo     
-          endif   !unres_pdb
-         else     !itype(i+1).eq.ntyp1
-          if (unres_pdb) then
-C 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
-            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+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 !itype(i+1).eq.ntyp1
-        endif  !itype.eq.ntyp1
-      enddo
+   10 continue
+#ifdef DEBUG
+      write (iout,'(a,i5)') ' Number of residues found: ',ires
+#endif
+      if (ires.eq.0) return
 C Calculate the CM of the last side chain.
+      if (iii.gt.0)  then
       if (unres_pdb) then
         do j=1,3
           dc(j,ires)=sccor(j,iii)
@@ -788,31 +430,19 @@ C Calculate the CM of the last side chain.
       else
         call sccenter(ires,iii,sccor)
       endif
+      endif
+      nres=ires
       nsup=nres
       nstart_sup=1
       if (itype(nres).ne.10) then
         nres=nres+1
         itype(nres)=ntyp1
-        if (unres_pdb) then
-C 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
         do j=1,3
-          dcj=(c(j,nres-2)-c(j,nres-3))/2.0
-            if (dcj.eq.0) dcj=1.23591524223
+          dcj=c(j,nres-2)-c(j,nres-3)
           c(j,nres)=c(j,nres-1)+dcj
           c(j,2*nres)=c(j,nres)
         enddo
       endif
-      endif
       do i=2,nres-1
         do j=1,3
           c(j,i+nres)=dc(j,i)
@@ -834,24 +464,18 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue
             e2(3)=0.0d0
           endif
           do j=1,3
-            c(j,1)=c(j,2)-1.9d0*e2(j)
+            c(j,1)=c(j,2)-3.8d0*e2(j)
           enddo
         else
         do j=1,3
-          dcj=(c(j,4)-c(j,3))/2.0
+          dcj=c(j,4)-c(j,3)
           c(j,1)=c(j,2)-dcj
           c(j,nres+1)=c(j,1)
         enddo
         endif
       endif
-C Copy the coordinates to reference coordinates
-c      do i=1,2*nres
-c        do j=1,3
-c          cref(j,i)=c(j,i)
-c        enddo
-c      enddo
 C Calculate internal coordinates.
-      if (out_template_coord) then
+      if (lprn) then
       write (iout,'(/a)') 
      &  "Cartesian coordinates of the reference structure"
       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') 
@@ -863,7 +487,16 @@ C Calculate internal coordinates.
       enddo
       endif
 C Calculate internal coordinates.
-      call int_from_cart(.true.,.true.)
+      if(me.eq.king.or..not.out1file)then
+       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)') 
+     &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
+     &    (c(j,nres+ires),j=1,3)
+       enddo
+      endif
+      call int_from_cart(.true.,.false.)
       call sc_loc_geom(.false.)
       do i=1,nres
         thetaref(i)=theta(i)
@@ -883,18 +516,16 @@ C Calculate internal coordinates.
 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 c     &   vbld_inv(i+nres)
       enddo
-      do i=1,nres
-        do j=1,3
-          cref(j,i)=c(j,i)
-          cref(j,i+nres)=c(j,i+nres)
-        enddo
-      enddo
+c      call chainbuild
+C Copy the coordinates to reference coordinates
       do i=1,2*nres
         do j=1,3
+          cref(j,i)=c(j,i)
           chomo(j,i,k)=c(j,i)
         enddo
       enddo
 
+
+      ishift_pdb=ishift
       return
       end
-      
index 8f4d05d..0320484 100644 (file)
@@ -177,8 +177,8 @@ c      call readi(controlcard,'IZ_SC',iz_sc,0)
       dccart=(index(controlcard,'CART').gt.0)
       overlapsc=(index(controlcard,'OVERLAP').gt.0)
       overlapsc=.not.overlapsc
-      searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
-      searchsc=.not.searchsc
+      searchsc=(index(controlcard,'SEARCHSC').gt.0 .and.
+     &  index(controlcard,'NOSEARCHSC').eq.0)
       sideadd=(index(controlcard,'SIDEADD').gt.0)
       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
       mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
index b467d9c..dc976b8 100644 (file)
@@ -19,6 +19,7 @@ C Loop over chain permutations
       if (iz_sc.lt.2) then
         do ichain=1,nchain
           indchain=tabpermchain(ichain,iperm)
+#define DEBUG
 #ifdef DEBUG
           write (iout,*) "ichain",ichain," indchain",indchain
           write (iout,*) "chain_border",chain_border(1,ichain),
@@ -199,8 +200,9 @@ c------------------------------------------------------------------------
       enddo
       rmsside=dsqrt(rmsside/nnnn)
 #ifdef DEBUG
-      write (iout,*) iii,iref," nnnn",nnnn," rmsside",rmsside
+      write (iout,*) "nnnn",nnnn," rmsside",rmsside
 #endif
+#undef DEBUG
       rmscalc_side=rmsside
       return
       end
index 75b7211..201e20f 100644 (file)
@@ -156,6 +156,7 @@ c     Put the original weights back to calculate the full energy
       mask_side=1
 crc      n_fun=n_fun+1
 ct      write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime
+      call chainbuild_extconf
       return
       end
 
index 09b6877..c8a3a0d 100644 (file)
@@ -525,7 +525,7 @@ c      save licznik
       integer IERROR
       integer i,j,k,l,ind,ind1,m,ii,iti,it,nzero,innt,inct
       integer ichain,nind
-      logical lprn /.true./
+      logical lprn /.false./
       double precision dtdi,gamvec(MAXRES2)
       common /syfek/ gamvec
 #ifndef FIVEDIAG
index f556eb6..08caade 100644 (file)
@@ -283,9 +283,9 @@ crc overlap test
         endif
 
         if (overlapsc) then 
-c          print *, 'Calling OVERLAP_SC'
+          write (iout,*) 'Calling OVERLAP_SC'
           call overlap_sc(fail)
-c          print *,"After overlap_sc"
+          write (iout,*) "After overlap_sc"
         endif 
 
         if (searchsc) then 
index e667382..5280f0a 100644 (file)
@@ -1,9 +1,9 @@
 BIN = ~/bin
 FC = ftn
-OPT = -mcmodel=medium -shared-intel -O3 -dynamic
+#OPT = -mcmodel=medium -shared-intel -O3 -dynamic
 #OPT = -O3 -intel-static -mcmodel=medium 
 #OPT = -O3 -ip -w 
-#OPT = -g -CB -mcmodel=medium -shared-intel -dynamic
+OPT = -g -CB -mcmodel=medium -shared-intel -dynamic
 FFLAGS = ${OPT} -c -I. -I./include_unres -I$(INSTALL_DIR)/include
 LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdrf.a
 
index b8ce4f4..cde2738 100644 (file)
@@ -1,10 +1,11 @@
       subroutine readpdb
 C Read the PDB file and convert the peptide geometry into virtual-chain 
 C geometry.
-      implicit none
+      implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'DIMENSIONS.ZSCOPT'
-      include 'COMMON.CONTROL'
+      include 'DIMENSIONS.FREE'
+      include 'COMMON.FRAG'
       include 'COMMON.LOCAL'
       include 'COMMON.VAR'
       include 'COMMON.CHAIN'
@@ -12,146 +13,176 @@ C geometry.
       include 'COMMON.IOUNITS'
       include 'COMMON.GEO'
       include 'COMMON.NAMES'
-      include 'COMMON.SBRIDGE'
-      character*3 seq,atom,res
+      include 'COMMON.CONTROL'
+      integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity
+      logical lprn /.false./,fail
+      double precision e1(3),e2(3),e3(3)
+      double precision dcj,efree_temp
+      character*3 seq,res
+      character*5 atom
       character*80 card
-      double precision sccor(3,50)
-      integer i,j,iii,ibeg,ishift,ishift1,ity,ires,ires_old
-      double precision dcj
-      integer rescode,kkk,lll,icha,cou,kupa,iprzes
+      double precision sccor(3,20)
+      integer rescode
+      efree_temp=0.0d0
       ibeg=1
       ishift1=0
+      ishift=0
+c      write (2,*) "UNRES_PDB",unres_pdb
+      ires=0
+      ires_old=0
+      iii=0
+      lsecondary=.false.
+      nhfrag=0
+      nbfrag=0
       do
         read (ipdbin,'(a80)',end=10) card
-        if (card(:3).eq.'END') then
-          goto 10
-        else if (card(:3).eq.'TER') then
-C End current chain
-c          ires_old=ires+1 
-          ires_old=ires+2
-          itype(ires_old-1)=ntyp1 
-          itype(ires_old)=ntyp1
-          ibeg=2
-c          write (iout,*) "Chain ended",ires,ishift,ires_old
-          call sccenter(ires,iii,sccor)
+c        write (iout,'(a)') card
+        if (card(:5).eq.'HELIX') then
+         nhfrag=nhfrag+1
+         lsecondary=.true.
+         read(card(22:25),*) hfrag(1,nhfrag)
+         read(card(34:37),*) hfrag(2,nhfrag)
+        endif
+        if (card(:5).eq.'SHEET') then
+         nbfrag=nbfrag+1
+         lsecondary=.true.
+         read(card(24:26),*) bfrag(1,nbfrag)
+         read(card(35:37),*) bfrag(2,nbfrag)
+crc----------------------------------------
+crc  to be corrected !!!
+         bfrag(3,nbfrag)=bfrag(1,nbfrag)
+         bfrag(4,nbfrag)=bfrag(2,nbfrag)
+crc----------------------------------------
         endif
+        if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
+c Read free energy
+        if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
 C Fish out the ATOM cards.
         if (index(card(1:4),'ATOM').gt.0) then  
-          read (card(14:16),'(a3)') atom
-          if (atom.eq.'CA' .or. atom.eq.'CH3') then
+          read (card(12:16),*) atom
+c          write (iout,*) "! ",atom," !",ires
+c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
+          read (card(23:26),*) ires
+          read (card(18:20),'(a3)') res
+c          write (iout,*) "ires",ires,ires-ishift+ishift1,
+c     &      " ires_old",ires_old
+c          write (iout,*) "ishift",ishift," ishift1",ishift1
+c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
+          if (ires-ishift+ishift1.ne.ires_old) then
 C Calculate the CM of the preceding residue.
+c            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
             if (ibeg.eq.0) then
-              call sccenter(ires,iii,sccor)
+c              write (iout,*) "Calculating sidechain center iii",iii
+              if (unres_pdb) then
+                do j=1,3
+                  dc(j,ires)=sccor(j,iii)
+                enddo
+              else
+                call sccenter(ires_old,iii,sccor)
+              endif
+              iii=0
             endif
 C Start new residue.
-c            write (iout,'(a80)') card
-            read (card(23:26),*) ires
-            read (card(18:20),'(a3)') res
-            if (ibeg.eq.1) then
+            if (res.eq.'Cl-' .or. res.eq.'Na+') then
+              ires=ires_old
+              cycle
+            else if (ibeg.eq.1) then
+c              write (iout,*) "BEG ires",ires
               ishift=ires-1
               if (res.ne.'GLY' .and. res.ne. 'ACE') then
                 ishift=ishift-1
                 itype(1)=ntyp1
               endif
-c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
+              ires=ires-ishift+ishift1
+              ires_old=ires
+c              write (iout,*) "ishift",ishift," ires",ires,
+c     &         " ires_old",ires_old
               ibeg=0          
-            else if (ibeg.eq.2) then
-c Start a new chain
-              ishift=-ires_old+ires-1
-c              write (iout,*) "New chain started",ires,ishift
-              ibeg=0
+            else
+              ishift=ishift-(ires-ishift+ishift1-ires_old-1)
+              ires=ires-ishift+ishift1
+              ires_old=ires
             endif
-            ires=ires-ishift
-c            write (2,*) "ires",ires," ishift",ishift
-            if (res.eq.'ACE') then
-              ity=10
+            if (res.eq.'ACE' .or. res.eq.'NHE') then
+              itype(ires)=10
             else
               itype(ires)=rescode(ires,res,0)
             endif
+          else
+            ires=ires-ishift+ishift1
+          endif
+c          write (iout,*) "ires_old",ires_old," ires",ires
+          if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+c            ishift1=ishift1+1
+          endif
+c          write (2,*) "ires",ires," res ",res," ity",ity
+          if (atom.eq.'CA' .or. atom.eq.'CH3' .or. 
+     &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
-            read(card(61:66),*) bfac(ires)
-c            write (iout,'(2i3,2x,a,3f8.3,5x,f8.3)') 
-c     &       ires,itype(ires),res,(c(j,ires),j=1,3),bfac(ires)
-            iii=1
+c            write (iout,*) "backbone ",atom 
+#ifdef DEBUG
+            write (iout,'(2i3,2x,a,3f8.3)') 
+     &      ires,itype(ires),res,(c(j,ires),j=1,3)
+#endif
+            iii=iii+1
             do j=1,3
               sccor(j,iii)=c(j,ires)
             enddo
-          else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and.
-     &             atom(1:1).ne.'Q' .and. atom(1:2).ne.'1H' .and.
-     &             atom(1:2).ne.'2H' .and. atom(1:2).ne.'3H' .and.
-     &             atom.ne.'N  ' .and. atom.ne.'C   ' .and.
-     &             atom.ne.'OXT' ) then
+            if (ishift.ne.0) then
+              ires_ca=ires+ishift-ishift1
+            else
+              ires_ca=ires
+            endif
+c            write (*,*) card(23:27),ires,itype(ires)
+          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.
+     &             atom.ne.'OXT' .and. atom(:2).ne.'3H') then
+c            write (iout,*) "sidechain ",atom
             iii=iii+1
-c            write (iout,*) res,ires,iii,atom
             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
-c            write (iout,'(3f8.3)') (sccor(j,iii),j=1,3)
           endif
         endif
       enddo
-   10 write (iout,'(a,i5)') ' Nres: ',ires
-C Calculate dummy residue coordinates inside the "chain" of a multichain
-C system
-      nres=ires
-      do i=2,nres-1
-c        write (iout,*) i,itype(i)
-
-        if (itype(i).eq.ntyp1) then
-         if (itype(i+1).eq.ntyp1) then
-C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
-C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
-C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
-C           if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the last dummy residue
-C            call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
-C            if (fail) then
-C              e2(1)=0.0d0
-C              e2(2)=1.0d0
-C              e2(3)=0.0d0
-C            endif !fail
-C            do j=1,3
-C             c(j,i)=c(j,i-1)-1.9d0*e2(j)
-C            enddo
-C           else   !unres_pdb
-           do j=1,3
-             dcj=(c(j,i-2)-c(j,i-3))/2.0
-             c(j,i)=c(j,i-1)+dcj
-             c(j,nres+i)=c(j,i)
-           enddo     
-C          endif   !unres_pdb
-         else     !itype(i+1).eq.ntyp1
-C          if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the first dummy residue
-C            call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
-C            if (fail) then
-C              e2(1)=0.0d0
-C              e2(2)=1.0d0
-C              e2(3)=0.0d0
-C            endif
-C            do j=1,3
-C              c(j,i)=c(j,i+1)-1.9d0*e2(j)
-C            enddo
-C          else !unres_pdb
-           do j=1,3
-            dcj=(c(j,i+3)-c(j,i+2))/2.0
-            c(j,i)=c(j,i+1)-dcj
-            c(j,nres+i)=c(j,i)
-           enddo
-C          endif !unres_pdb
-         endif !itype(i+1).eq.ntyp1
-        endif  !itype.eq.ntyp1
-      enddo
+   10 continue
+#ifdef DEBUG
+      write (iout,'(a,i5)') ' Number of residues found: ',ires
+#endif
+      if (ires.eq.0) return
 C Calculate the CM of the last side chain.
-      call sccenter(ires,iii,sccor)
+      if (iii.gt.0)  then
+      if (unres_pdb) then
+        do j=1,3
+          dc(j,ires)=sccor(j,iii)
+        enddo
+      else
+        call sccenter(ires,iii,sccor)
+      endif
+      endif
+      nres=ires
       nsup=nres
       nstart_sup=1
       if (itype(nres).ne.10) then
         nres=nres+1
         itype(nres)=ntyp1
+        if (unres_pdb) then
+C 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)-3.8d0*e2(j)
+          enddo
+        else
         do j=1,3
-          dcj=(c(j,nres-2)-c(j,nres-3))/2.0
+          dcj=c(j,nres-2)-c(j,nres-3)
           c(j,nres)=c(j,nres-1)+dcj
           c(j,2*nres)=c(j,nres)
         enddo
+        endif
       endif
       do i=2,nres-1
         do j=1,3
@@ -165,21 +196,57 @@ C Calculate the CM of the last side chain.
       if (itype(1).eq.ntyp1) then
         nsup=nsup-1
         nstart_sup=2
+        if (unres_pdb) then
+C 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)-3.8d0*e2(j)
+          enddo
+        else
         do j=1,3
-          dcj=(c(j,4)-c(j,3))/2.0
+          dcj=c(j,4)-c(j,3)
           c(j,1)=c(j,2)-dcj
           c(j,nres+1)=c(j,1)
         enddo
+        endif
       endif
+C Copy the coordinates to reference coordinates
+c      do i=1,2*nres
+c        do j=1,3
+c          cref(j,i)=c(j,i)
+c        enddo
+c      enddo
 C Calculate internal coordinates.
-      write (iout,100)
+      if (lprn) then
+      write (iout,'(/a)') 
+     &  "Cartesian coordinates of the reference structure"
+      write (iout,'(a,3(3x,a5),5x,3(3x,a5))') 
+     & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
       do ires=1,nres
+        write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') 
+     &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
+     &    (c(j,ires+nres),j=1,3)
+      enddo
+      endif
+C Calculate internal coordinates.
+       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)') 
      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
      &    (c(j,nres+ires),j=1,3)
-      enddo
+       enddo
       call int_from_cart(.true.,.false.)
-      call flush(iout)
+      call sc_loc_geom(.false.)
+      do i=1,nres
+        thetaref(i)=theta(i)
+        phiref(i)=phi(i)
+      enddo
       do i=1,nres-1
         do j=1,3
           dc(j,i)=c(j,i+1)-c(j,i)
@@ -196,24 +263,30 @@ c     &   vbld_inv(i+nres)
       enddo
 c      call chainbuild
 C Copy the coordinates to reference coordinates
-      do i=1,nres
+      do i=1,2*nres
         do j=1,3
           cref(j,i)=c(j,i)
-          cref(j,i+nres)=c(j,i+nres)
         enddo
       enddo
-  100 format ('Residue    alpha-carbon coordinates    ',
-     &          '     centroid coordinates'/
-     1          '         ', 6X,'X',7X,'Y',7X,'Z',
-     &                          12X,'X',7X,'Y',7X,'Z')
-  110 format (a,'(',i3,')',6f12.5)
 
+
+      do j=1,nbfrag     
+        do i=1,4                                                       
+         bfrag(i,j)=bfrag(i,j)-ishift
+        enddo
+      enddo
+
+      do j=1,nhfrag
+        do i=1,2
+         hfrag(i,j)=hfrag(i,j)-ishift
+        enddo
+      enddo
       ishift_pdb=ishift
       return
       end
 c---------------------------------------------------------------------------
       subroutine int_from_cart(lside,lprn)
-      implicit none
+      implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'DIMENSIONS.ZSCOPT'
       include 'COMMON.LOCAL'
@@ -223,56 +296,62 @@ c---------------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.GEO'
       include 'COMMON.NAMES'
-      character*3 seq,atom,res
+      include 'COMMON.CONTROL'
+      character*3 seq,res
+c      character*5 atom
       character*80 card
-      double precision sccor(3,50)
+      dimension sccor(3,20)
       integer rescode
-      double precision dist,alpha,beta,di
-      integer i,j,iti
       logical lside,lprn
-      if (lprn) then 
+       if (lprn) then 
         write (iout,'(/a)') 
      &  'Internal coordinates calculated from crystal structure.'
         if (lside) then 
           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
-     & '       Phi','    Dsc_id','       Dsc','     Alpha',
-     & '     Omega'
+     & '     Gamma','    Dsc_id','       Dsc','     Alpha',
+     & '     Beta '
         else 
           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
-     & '       Phi'
+     & '     Gamma'
         endif
-      endif
-      do i=2,nres
+       endif
+      do i=1,nres-1
         iti=itype(i)
-c        write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1)
-        if (itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1 .and.
-     &    (dist(i,i-1).lt.1.0D0 .or. dist(i,i-1).gt.6.0D0)) then
+        if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
-          stop
+ctest          stop
         endif
-        vbld(i)=dist(i-1,i)
-        vbld_inv(i)=1.0d0/vbld(i)
-        theta(i+1)=alpha(i-1,i,i+1)
+        vbld(i+1)=dist(i,i+1)
+        vbld_inv(i+1)=1.0d0/vbld(i+1)
+        if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
       enddo
-c      if (itype(1).eq.ntyp1) then
-c        do j=1,3
-c          c(j,1)=c(j,2)+(c(j,3)-c(j,4))
-c        enddo
-c      endif
-c      if (itype(nres).eq.ntyp1) then
-c        do j=1,3
-c          c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3))
-c        enddo
+c      if (unres_pdb) then
+c        if (itype(1).eq.ntyp1) then
+c          theta(3)=90.0d0*deg2rad
+c          phi(4)=180.0d0*deg2rad
+c          vbld(2)=3.8d0
+c          vbld_inv(2)=1.0d0/vbld(2)
+c        endif
+c        if (itype(nres).eq.ntyp1) then
+c          theta(nres)=90.0d0*deg2rad
+c          phi(nres)=180.0d0*deg2rad
+c          vbld(nres)=3.8d0
+c          vbld_inv(nres)=1.0d0/vbld(2)
+c        endif
 c      endif
       if (lside) then
         do i=2,nres-1
           do j=1,3
-            c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
+            c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
+     &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
           enddo
           iti=itype(i)
           di=dist(i,nres+i)
-           vbld(i+nres)=di
+C 10/03/12 Adam: Correction for zero SC-SC bond length
+          if (itype(i).ne.10 .and. itype(i).ne.ntyp1. and. di.eq.0.0d0)
+     &     di=dsc(itype(i))
+          vbld(i+nres)=di
           if (itype(i).ne.10) then
             vbld_inv(i+nres)=1.0d0/di
           else
@@ -282,41 +361,21 @@ c      endif
             alph(i)=alpha(nres+i,i,maxres2)
             omeg(i)=beta(nres+i,i,maxres2,i+1)
           endif
-          if (iti.ne.10) then
-            alph(i)=alpha(nres+i,i,maxres2)
-            omeg(i)=beta(nres+i,i,maxres2,i+1)
-          endif
-          if (lprn)
-     &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
-     &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,
-     &    rad2deg*alph(i),rad2deg*omeg(i)
+           if (lprn)
+     &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
+     &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
+     &     rad2deg*alph(i),rad2deg*omeg(i)
         enddo
       else if (lprn) then
         do i=2,nres
           iti=itype(i)
           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
-     &    rad2deg*theta(i),rad2deg*phi(i)
+     &     rad2deg*theta(i),rad2deg*phi(i)
         enddo
       endif
       return
       end
-c---------------------------------------------------------------------------
-      subroutine sccenter(ires,nscat,sccor)
-      implicit none
-      include 'DIMENSIONS'
-      include 'COMMON.CHAIN'
-      integer ires,nscat,i,j
-      double precision sccor(3,50),sccmj
-      do j=1,3
-        sccmj=0.0D0
-        do i=1,nscat
-          sccmj=sccmj+sccor(j,i) 
-        enddo
-        dc(j,ires)=sccmj/nscat
-      enddo
-      return
-      end
-c---------------------------------------------------------------------------
+c-------------------------------------------------------------------------------
       subroutine sc_loc_geom(lprn)
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
@@ -329,7 +388,6 @@ c---------------------------------------------------------------------------
       include 'COMMON.GEO'
       include 'COMMON.NAMES'
       include 'COMMON.CONTROL'
-      include 'COMMON.SETUP'
       double precision x_prime(3),y_prime(3),z_prime(3)
       logical lprn
       do i=1,nres-1
@@ -338,7 +396,7 @@ c---------------------------------------------------------------------------
         enddo
       enddo
       do i=2,nres-1
-        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+        if (itype(i).ne.10) then
           do j=1,3
             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
           enddo
@@ -358,7 +416,7 @@ c---------------------------------------------------------------------------
         sinfac2=0.5d0/(1.0d0-costtab(i+1))
         sinfac=dsqrt(sinfac2)
         it=itype(i)
-        if (it.ne.10 .and. itype(i).ne.ntyp1) then
+        if (it.ne.10) then
 c
 C  Compute the axes of tghe local cartesian coordinates system; store in
 c   x_prime, y_prime and z_prime 
@@ -372,10 +430,7 @@ c
           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
         enddo
-c        write (iout,*) "x_prime",(x_prime(j),j=1,3)
-c        write (iout,*) "y_prime",(y_prime(j),j=1,3)
         call vecpr(x_prime,y_prime,z_prime)
-c        write (iout,*) "z_prime",(z_prime(j),j=1,3)
 c
 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 C to local coordinate system. Store in xx, yy, zz.
@@ -399,16 +454,30 @@ c
         endif
       enddo
       if (lprn) then
-        write (iout,*) "xxref,yyref,zzref"
         do i=2,nres
           iti=itype(i)
-          write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),yyref(i),
-     &     zzref(i)
+          write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
+     &      yyref(i),zzref(i)
         enddo
       endif
       return
       end
 c---------------------------------------------------------------------------
+      subroutine sccenter(ires,nscat,sccor)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      dimension sccor(3,20)
+      do j=1,3
+        sccmj=0.0D0
+        do i=1,nscat
+          sccmj=sccmj+sccor(j,i) 
+        enddo
+        dc(j,ires)=sccmj/nscat
+      enddo
+      return
+      end
+c---------------------------------------------------------------------------
       subroutine bond_regular
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'   
@@ -420,19 +489,19 @@ c---------------------------------------------------------------------------
       do i=1,nres-1
        vbld(i+1)=vbl
        vbld_inv(i+1)=1.0d0/vbld(i+1)
-       vbld(i+1+nres)=dsc(iabs(itype(i+1)))
-       vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1)))
+       vbld(i+1+nres)=dsc(itype(i+1))
+       vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
 c       print *,vbld(i+1),vbld(i+1+nres)
       enddo
       return
       end
-c---------------------------------------------------------------------------
       subroutine readpdb_template(k)
-C Read the PDB file for read_constr_homology with read2sigma
+C Read the PDB file with gaps for read_constr_homology with read2sigma
 C and convert the peptide geometry into virtual-chain geometry.
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.FREE'
       include 'COMMON.LOCAL'
       include 'COMMON.VAR'
       include 'COMMON.CHAIN'
@@ -441,7 +510,6 @@ C and convert the peptide geometry into virtual-chain geometry.
       include 'COMMON.GEO'
       include 'COMMON.NAMES'
       include 'COMMON.CONTROL'
-      include 'COMMON.SETUP'
       integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity
       logical lprn /.false./,fail
       double precision e1(3),e2(3),e3(3)
@@ -450,10 +518,8 @@ C and convert the peptide geometry into virtual-chain geometry.
       character*5 atom
       character*80 card
       double precision sccor(3,20)
-      integer rescode,iterter(maxres)
-      do i=1,maxres
-         iterter(i)=0
-      enddo
+      integer rescode
+      efree_temp=0.0d0
       ibeg=1
       ishift1=0
       ishift=0
@@ -464,27 +530,10 @@ c      write (2,*) "UNRES_PDB",unres_pdb
       lsecondary=.false.
       nhfrag=0
       nbfrag=0
-      do
+      do 
         read (ipdbin,'(a80)',end=10) card
-        if (card(:3).eq.'END') then
-          goto 10
-        else if (card(:3).eq.'TER') then
-C End current chain
-          ires_old=ires+2
-          itype(ires_old-1)=ntyp1 
-          iterter(ires_old-1)=1
-          itype(ires_old)=ntyp1
-          iterter(ires_old)=1
-          ibeg=2
-c          write (iout,*) "Chain ended",ires,ishift,ires_old
-          if (unres_pdb) then
-            do j=1,3
-              dc(j,ires)=sccor(j,iii)
-            enddo
-          else 
-            call sccenter(ires,iii,sccor)
-          endif
-        endif
+c        write (iout,'(a)') card
+        if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
 C Fish out the ATOM cards.
         if (index(card(1:4),'ATOM').gt.0) then  
           read (card(12:16),*) atom
@@ -498,7 +547,9 @@ c          write (iout,*) "ishift",ishift," ishift1",ishift1
 c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
           if (ires-ishift+ishift1.ne.ires_old) then
 C Calculate the CM of the preceding residue.
+c            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
             if (ibeg.eq.0) then
+c              write (iout,*) "Calculating sidechain center iii",iii
               if (unres_pdb) then
                 do j=1,3
                   dc(j,ires)=sccor(j,iii)
@@ -523,13 +574,6 @@ c              write (iout,*) "BEG ires",ires
               ires_old=ires
 c              write (iout,*) "ishift",ishift," ires",ires,
 c     &         " ires_old",ires_old
-c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
-              ibeg=0          
-            else if (ibeg.eq.2) then
-c Start a new chain
-              ishift=-ires_old+ires-1
-              ires=ires_old+1
-c              write (iout,*) "New chain started",ires,ishift
               ibeg=0          
             else
               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
@@ -545,14 +589,14 @@ c              write (iout,*) "New chain started",ires,ishift
             ires=ires-ishift+ishift1
           endif
 c          write (iout,*) "ires_old",ires_old," ires",ires
-c          if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+          if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 c            ishift1=ishift1+1
-c          endif
+          endif
 c          write (2,*) "ires",ires," res ",res," ity",ity
           if (atom.eq.'CA' .or. atom.eq.'CH3' .or. 
      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
-c            write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
+c            write (iout,*) "backbone ",atom 
 #ifdef DEBUG
             write (iout,'(2i3,2x,a,3f8.3)') 
      &      ires,itype(ires),res,(c(j,ires),j=1,3)
@@ -577,60 +621,13 @@ c            write (iout,*) "sidechain ",atom
           endif
         endif
       enddo
-   10 write (iout,'(a,i5)') ' Nres: ',ires
-C Calculate dummy residue coordinates inside the "chain" of a multichain
-C system
-      nres=ires
-      do i=2,nres-1
-c        write (iout,*) i,itype(i),itype(i+1)
-        if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
-         if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
-C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
-C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
-C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
-           if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the last dummy residue
-            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
-            do j=1,3
-             c(j,i)=c(j,i-1)-1.9d0*e2(j)
-            enddo
-           else   !unres_pdb
-           do j=1,3
-             dcj=(c(j,i-2)-c(j,i-3))/2.0
-            if (dcj.eq.0) dcj=1.23591524223
-             c(j,i)=c(j,i-1)+dcj
-             c(j,nres+i)=c(j,i)
-           enddo     
-          endif   !unres_pdb
-         else     !itype(i+1).eq.ntyp1
-          if (unres_pdb) then
-C 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
-            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+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 !itype(i+1).eq.ntyp1
-        endif  !itype.eq.ntyp1
-      enddo
+   10 continue
+#ifdef DEBUG
+      write (iout,'(a,i5)') ' Number of residues found: ',ires
+#endif
+      if (ires.eq.0) return
 C Calculate the CM of the last side chain.
+      if (iii.gt.0)  then
       if (unres_pdb) then
         do j=1,3
           dc(j,ires)=sccor(j,iii)
@@ -638,31 +635,19 @@ C Calculate the CM of the last side chain.
       else
         call sccenter(ires,iii,sccor)
       endif
+      endif
+      nres=ires
       nsup=nres
       nstart_sup=1
       if (itype(nres).ne.10) then
         nres=nres+1
         itype(nres)=ntyp1
-        if (unres_pdb) then
-C 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
         do j=1,3
-          dcj=(c(j,nres-2)-c(j,nres-3))/2.0
-            if (dcj.eq.0) dcj=1.23591524223
+          dcj=c(j,nres-2)-c(j,nres-3)
           c(j,nres)=c(j,nres-1)+dcj
           c(j,2*nres)=c(j,nres)
         enddo
       endif
-      endif
       do i=2,nres-1
         do j=1,3
           c(j,i+nres)=dc(j,i)
@@ -684,24 +669,18 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue
             e2(3)=0.0d0
           endif
           do j=1,3
-            c(j,1)=c(j,2)-1.9d0*e2(j)
+            c(j,1)=c(j,2)-3.8d0*e2(j)
           enddo
         else
         do j=1,3
-          dcj=(c(j,4)-c(j,3))/2.0
+          dcj=c(j,4)-c(j,3)
           c(j,1)=c(j,2)-dcj
           c(j,nres+1)=c(j,1)
         enddo
         endif
       endif
-C Copy the coordinates to reference coordinates
-c      do i=1,2*nres
-c        do j=1,3
-c          cref(j,i)=c(j,i)
-c        enddo
-c      enddo
 C Calculate internal coordinates.
-      if (out_template_coord) then
+      if (lprn) then
       write (iout,'(/a)') 
      &  "Cartesian coordinates of the reference structure"
       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') 
@@ -713,9 +692,15 @@ C Calculate internal coordinates.
       enddo
       endif
 C Calculate internal coordinates.
-c      call int_from_cart1(.false.)
-      call int_from_cart(.true.,.true.)
-      call sc_loc_geom(.true.)
+       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)') 
+     &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
+     &    (c(j,nres+ires),j=1,3)
+       enddo
+      call int_from_cart(.true.,.false.)
+      call sc_loc_geom(.false.)
       do i=1,nres
         thetaref(i)=theta(i)
         phiref(i)=phi(i)
@@ -734,19 +719,16 @@ c      call int_from_cart1(.false.)
 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 c     &   vbld_inv(i+nres)
       enddo
-      do i=1,nres
-        do j=1,3
-          cref(j,i)=c(j,i)
-          cref(j,i+nres)=c(j,i+nres)
-        enddo
-      enddo
+c      call chainbuild
+C Copy the coordinates to reference coordinates
       do i=1,2*nres
         do j=1,3
+          cref(j,i)=c(j,i)
           chomo(j,i,k)=c(j,i)
         enddo
       enddo
 
+
+      ishift_pdb=ishift
       return
       end
-      
-