corrrection in energies (no trash)

diff --git a/source/cluster/CMakeLists.txt b/source/cluster/CMakeLists.txt
index 1f39d0f..9f80cd2 100644 (file)
--- a/source/cluster/CMakeLists.txt
+++ b/source/cluster/CMakeLists.txt
@@ -47,7 +47,7 @@ set(UNRES_CLUSTER_WHAM_SRC0
 if (Fortran_COMPILER_NAME STREQUAL "ifort")
   set (CMAKE_Fortran_FLAGS_RELEASE " ")
   set (CMAKE_Fortran_FLAGS_DEBUG   "-O0 -g ")     
-  set(FFLAGS0 "-fpp -mcmodel=medium -shared-intel -ip " )
+  set(FFLAGS0 "-CB -g -fpp -mcmodel=medium -shared-intel -ip " )
 elseif (Fortran_COMPILER_NAME STREQUAL "gfortran")
   set(FFLAGS0 "-std=legacy -mcmodel=medium " )
 elseif (Fortran_COMPILER_NAME STREQUAL "pgf90")

diff --git a/source/cluster/io_clust.f90 b/source/cluster/io_clust.f90
index db74bf6..7570389 100644 (file)
--- a/source/cluster/io_clust.f90
+++ b/source/cluster/io_clust.f90
@@ -4,7 +4,7 @@
       use io_units
 !      use names
       use io_base !, only: ilen
-      use geometry_data, only: nres,c
+      use geometry_data, only: nres,c,nres_molec,boxxsize,boxysize,boxzsize
       use energy_data, only: nnt,nct,nss
       use control_data, only: lside
       implicit none
@@ -865,7 +865,7 @@
       subroutine add_new_cconf(jjj,jj,jj_old,icount,Next)
 
       use geometry_data, only: vbld,rad2deg,theta,phi,alph,omeg,deg2rad
-      use energy_data, only: itel,itype,dsc,max_ene
+      use energy_data, only: itel,itype,dsc,max_ene,molnum
       use control_data, only: symetr
       use geometry, only: int_from_cart1
 !      implicit none
@@ -882,13 +882,18 @@
 !      include "COMMON.SBRIDGE"
 !      include "COMMON.GEO"
       integer :: i,j,jj,jjj,jj_old,icount,k,kk,l,ii,ib,&
-        nn,nn1,inan,Next,itj,chalen
+        nn,nn1,inan,Next,itj,chalen,mnum
       real(kind=8) :: etot,energia(0:max_ene)
       jjj=jjj+1
       chalen=int((nct-nnt+2)/symetr)
       call int_from_cart1(.false.)
       do j=nnt+1,nct
-        if (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0) then
+        mnum=molnum(j)
+        if (mnum.eq.5) cycle
+        if (itype(j,mnum).eq.ntyp1_molec(mnum)) cycle
+        if (itype(j+1,molnum(j+1)).eq.ntyp1_molec(molnum(j+1))) cycle
+
+        if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0).and.(mnum.ne.5)) then
          if (j.gt.2) then
           if (itel(j).ne.0 .and. itel(j-1).ne.0) then
           write (iout,*) "Conformation",jjj,jj+1
@@ -912,8 +917,11 @@
         endif
       enddo
       do j=nnt,nct
-        itj=itype(j)
-        if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(iabs(itj))) &
+        mnum=molnum(j)
+        itj=itype(j,mnum)
+        if (mnum.eq.5) cycle
+        if (itype(j,1).ne.10 .and. (itype(j,mnum).ne.ntyp1_molec(mnum)) & 
+            .and.  (vbld(nres+j)-dsc(iabs(itj))) &
                                   .gt.2.0d0) then
           write (iout,*) "Conformation",jjj,jj+1
           write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j)
@@ -1248,7 +1256,8 @@
       use energy_data, only: rescale_mode,distchainmax,ipot !,temp0
       use control_data, only: titel,outpdb,outmol2,refstr,pdbref,&
                  iscode,symetr,punch_dist,print_dist,nstart,nend,&
-                 caonly,iopt,efree,lprint_cart,lprint_int
+                 caonly,iopt,efree,lprint_cart,lprint_int,rlamb_ele,&
+                 r_cut_ele
 !      implicit none
 !      include 'DIMENSIONS'
 !      include 'sizesclu.dat'
@@ -1271,7 +1280,13 @@
       read (INP,'(a80)') titel
       call card_concat(controlcard,.true.)
 
-      call readi(controlcard,'NRES',nres,0)
+      call readi(controlcard,'NRES',nres_molec(1),0)
+      call readi(controlcard,"NRES_NUCL",nres_molec(2),0)
+      call readi(controlcard,"NRES_CAT",nres_molec(5),0)
+      nres=0
+      do i=1,5
+       nres=nres_molec(i)+nres
+      enddo
 
 !      call alloc_clust_arrays
       allocate(rcutoff(max_cut+1)) !(max_cut+1)
@@ -1296,6 +1311,13 @@
       call readi(controlcard,'SYM',symetr,1)
       write (iout,*) 'sym', symetr
       call readi(controlcard,'NSTART',nstart,0)
+      call reada(controlcard,'BOXX',boxxsize,100.0d0)
+      call reada(controlcard,'BOXY',boxysize,100.0d0)
+      call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+!      ions=index(controlcard,"IONS").gt.0
+      call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
+      call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
+      write(iout,*) "R_CUT_ELE=",r_cut_ele
       call readi(controlcard,'NEND',nend,0)
       call reada(controlcard,'ECUT',ecut,10.0d0)
       call reada(controlcard,'PROB',prob_limit,0.99d0)
@@ -1361,7 +1383,7 @@
       real(kind=8) :: x(6*nres) !(maxvar)
       integer  :: itype_pdb(nres) !(maxres)
 !      logical seq_comp
-      integer :: i,j,kkk
+      integer :: i,j,kkk,mnum
 !
 ! Body
 !
@@ -1372,7 +1394,8 @@
 !-----------------------------
       allocate(c(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
       allocate(dc(3,0:2*nres+2)) !(3,0:maxres2)
-      allocate(itype(nres+2)) !(maxres)
+      allocate(itype(nres+2,5)) !(maxres)
+      allocate(molnum(nres+1))
       allocate(itel(nres+2))
 
       do i=1,2*nres+2
@@ -1381,10 +1404,12 @@
           dc(j,i)=0
         enddo
       enddo
+      do j=1,5
       do i=1,nres+2
-        itype(i)=0
+        itype(i,j)=0
         itel(i)=0
       enddo
+      enddo
 !--------------------------
 ! Read weights of the subsequent energy terms.
       call card_concat(weightcard,.true.)
@@ -1411,6 +1436,8 @@
       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
+      call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
+      call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
       call reada(weightcard,'TEMP0',temp0,300.0d0)   !!! el
       if (index(weightcard,'SOFT').gt.0) ipot=6
 ! 12/1/95 Added weight for the multi-body term WCORR
@@ -1491,23 +1518,36 @@
         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
       endif
 ! Convert sequence to numeric code
-      do i=1,nres
-        itype(i)=rescode(i,sequence(i),iscode)
+      do i=1,nres_molec(1)
+        mnum=1
+        molnum(i)=1
+        itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
+      enddo
+      do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
+        mnum=2
+        molnum(i)=2
+        itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
+      enddo
+      do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
+        mnum=5
+        molnum(i)=5
+        itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
       enddo
       print *,nres
-      print '(20i4)',(itype(i),i=1,nres)
+      print '(20i4)',(itype(i,molnum(i)),i=1,nres)
 
-      do i=1,nres
+      do i=1,nres-1
 #ifdef PROCOR
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
+        if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))&
+            .or. itype(i+1,molnum(i+1)).eq.ntyp1_molec(molnum(i+1))) then
 #else
-        if (itype(i).eq.ntyp1) then
+        if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) then
 #endif
           itel(i)=0
 #ifdef PROCOR
-        else if (iabs(itype(i+1)).ne.20) then
+        else if (iabs(itype(i+1,1)).ne.20) then
 #else
-        else if (iabs(itype(i)).ne.20) then
+        else if (iabs(itype(i,1)).ne.20) then
 #endif
           itel(i)=1
         else
@@ -1516,7 +1556,7 @@
       enddo
       write (iout,*) "ITEL"
       do i=1,nres-1
-        write (iout,*) i,itype(i),itel(i)
+        write (iout,*) i,itype(i,molnum(i)),itel(i)
       enddo
 
       print *,'Call Read_Bridge.'
@@ -1524,8 +1564,8 @@
       nnt=1
       nct=nres
       print *,'NNT=',NNT,' NCT=',NCT
-      if (itype(1).eq.ntyp1) nnt=2
-      if (itype(nres).eq.ntyp1) nct=nct-1
+      if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
+      if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
       if (nstart.lt.nnt) nstart=nnt
       if (nend.gt.nct .or. nend.eq.0) nend=nct
       write (iout,*) "nstart",nstart," nend",nend
@@ -1681,6 +1721,33 @@ write(iout,*)"po setup var"
       open (isidep1,file=sidepname,status="old")
       call getenv('SCCORPAR',sccorname)
       open (isccor,file=sccorname,status="old")
+      call getenv('THETPAR_NUCL',thetname_nucl)
+      open (ithep_nucl,file=thetname_nucl,status='old')
+      call getenv('ROTPAR_NUCL',rotname_nucl)
+      open (irotam_nucl,file=rotname_nucl,status='old')
+      call getenv('TORPAR_NUCL',torname_nucl)
+      open (itorp_nucl,file=torname_nucl,status='old')
+      call getenv('TORDPAR_NUCL',tordname_nucl)
+      open (itordp_nucl,file=tordname_nucl,status='old')
+      call getenv('SIDEPAR_NUCL',sidename_nucl)
+      open (isidep_nucl,file=sidename_nucl,status='old')
+      call getenv('SIDEPAR_SCBASE',sidename_scbase)
+      open (isidep_scbase,file=sidename_scbase,status='old')
+      call getenv('PEPPAR_PEPBASE',pepname_pepbase)
+      open (isidep_pepbase,file=pepname_pepbase,status='old')
+      call getenv('SCPAR_PHOSPH',pepname_scpho)
+      open (isidep_scpho,file=pepname_scpho,status='old')
+      call getenv('PEPPAR_PHOSPH',pepname_peppho)
+      open (isidep_peppho,file=pepname_peppho,status='old')
+
+
+      call getenv('LIPTRANPAR',liptranname)
+      open (iliptranpar,file=liptranname,status='old')
+      call getenv('TUBEPAR',tubename)
+      open (itube,file=tubename,status='old')
+      call getenv('IONPAR',ionname)
+      open (iion,file=ionname,status='old')
+
 #ifndef OLDSCP
 !
 ! 8/9/01 In the newest version SCp interaction constants are read from a file
@@ -1696,7 +1763,7 @@ write(iout,*)"po setup var"
 !-----------------------------------------------------------------------------
       subroutine pdboutC(etot,rmsd,tytul)
 
-      use energy_data, only: ihpb,jhpb,itype
+      use energy_data, only: ihpb,jhpb,itype,molnum
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.CONTROL'
@@ -1712,16 +1779,28 @@ write(iout,*)"po setup var"
                                                'G','H','I','J'/)
       integer :: ica(nres)
       real(kind=8) :: etot,rmsd
-      integer :: iatom,ichain,ires,i,j,iti
-
+      integer :: iatom,ichain,ires,i,j,iti,mnum,mnum1,boxxxshift(3)
+      real(kind=8) :: boxxxx(3)
       write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
         ' ENERGY ',etot,' RMS ',rmsd
       iatom=0
       ichain=1
       ires=0
+      boxxxshift(1)=int(c(1,nnt)/boxxsize)
+!      if (c(1,nnt).lt.0) boxxxshift(1)=boxxxshift(1)-1
+      boxxxshift(2)=int(c(2,nnt)/boxzsize)
+!      if (c(2,nnt).lt.0) boxxxshift(2)=boxxxshift(2)-1
+      boxxxshift(3)=int(c(3,nnt)/boxzsize)
+!      if (c(3,nnt).lt.0) boxxxshift(3)=boxxxshift(3)-1
+
+      boxxxx(1)=boxxsize
+      boxxxx(2)=boxysize
+      boxxxx(3)=boxzsize
+
       do i=nnt,nct
-        iti=itype(i)
-        if (iti.eq.ntyp1) then
+        mnum=molnum(i)
+        iti=itype(i,mnum)
+        if (iti.eq.ntyp1_molec(mnum)) then
           ichain=ichain+1
           ires=0
           write (ipdb,'(a)') 'TER'
@@ -1729,27 +1808,40 @@ write(iout,*)"po setup var"
         ires=ires+1
         iatom=iatom+1
         ica(i)=iatom
-        write (ipdb,10) iatom,restyp(iti),chainid(ichain),&
+        if (mnum.eq.5) then
+           do j=1,3
+            if ((c(j,i).gt.0).and.(c(j,nnt).lt.0)) then
+           c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)-1)*boxxxx(j)
+            else if ((c(j,i).lt.0).and.(c(j,nnt).gt.0)) then
+           c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)+1)*boxxxx(j)
+             else
+           c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j))*boxxxx(j)
+            endif
+           enddo
+        endif
+        write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
            ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
-        if (iti.ne.10) then
+        if ((iti.ne.10).and.(mnum.ne.5)) then
           iatom=iatom+1
-          write (ipdb,20) iatom,restyp(iti),chainid(ichain),&
+          write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
             ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i)
         endif
         endif
       enddo
       write (ipdb,'(a)') 'TER'
       do i=nnt,nct-1
-        if (itype(i).eq.ntyp1) cycle
-        if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
+        mnum=molnum(i)
+        mnum1=molnum(i+1)
+        if ((itype(i,mnum).eq.ntyp1).or.(mnum.eq.5)) cycle
+        if (itype(i,mnum).eq.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
           write (ipdb,30) ica(i),ica(i+1)
-        else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
+        else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
           write (ipdb,30) ica(i),ica(i+1),ica(i)+1
-        else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
+        else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) then
           write (ipdb,30) ica(i),ica(i)+1
         endif
       enddo
-      if (itype(nct).ne.10) then
+      if ((itype(nct,molnum(nct)).ne.10).and.(molnum(nct).ne.5)) then
         write (ipdb,30) ica(nct),ica(nct)+1
       endif
       do i=1,nss

diff --git a/source/unres/io.f90 b/source/unres/io.f90
index 05f3585..ab5dc2a 100644 (file)
--- a/source/unres/io.f90
+++ b/source/unres/io.f90
@@ -769,19 +769,19 @@
       call reada(weightcard,'WTURN6',wturn6,1.0D0)
       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
-      call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,1.0D0)
-      call reada(weightcard,'WELPP',welpp,1.0d0)
-      call reada(weightcard,'WVDWPSB',wvdwpsb,1.0d0)
-      call reada(weightcard,'WELPSB',welpsb,1.0D0)
-      call reada(weightcard,'WVDWSB',wvdwsb,1.0d0)
-      call reada(weightcard,'WELSB',welsb,1.0D0)
-      call reada(weightcard,'WBOND_NUCL',wbond_nucl,1.0D0)
-      call reada(weightcard,'WANG_NUCL',wang_nucl,1.0D0)
-      call reada(weightcard,'WSBLOC',wsbloc,1.0D0)
-      call reada(weightcard,'WTOR_NUCL',wtor_nucl,1.0D0)
-      call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,1.0D0)
-      call reada(weightcard,'WCORR_NUCL',wcorr_nucl,1.0D0)
-      call reada(weightcard,'WCORR3_NUCL',wcorr3_nucl,1.0D0)
+      call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,0.0D0)
+      call reada(weightcard,'WELPP',welpp,0.0d0)
+      call reada(weightcard,'WVDWPSB',wvdwpsb,0.0d0)
+      call reada(weightcard,'WELPSB',welpsb,0.0D0)
+      call reada(weightcard,'WVDWSB',wvdwsb,0.0d0)
+      call reada(weightcard,'WELSB',welsb,0.0D0)
+      call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
+      call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0)
+      call reada(weightcard,'WSBLOC',wsbloc,0.0D0)
+      call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0)
+      call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0)
+      call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0)
+      call reada(weightcard,'WCORR3_NUCL',wcorr3_nucl,0.0D0)
       call reada(weightcard,'WBOND',wbond,1.0D0)
       call reada(weightcard,'WTOR',wtor,1.0D0)
       call reada(weightcard,'WTORD',wtor_d,1.0D0)
@@ -796,14 +796,14 @@
       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
       call reada(weightcard,'TEMP0',temp0,300.0d0)
-      call reada(weightcard,'WSCBASE',wscbase,1.0D0)
+      call reada(weightcard,'WSCBASE',wscbase,0.0D0)
       if (index(weightcard,'SOFT').gt.0) ipot=6
-      call reada(weightcard,'WBOND_NUCL',wbond_nucl,1.0D0)
+      call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
       call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
       call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
       call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
-      call reada(weightcard,'WSCPHO',wscpho,1.0d0)
-      call reada(weightcard,'WPEPPHO',wpeppho,1.0d0)
+      call reada(weightcard,'WSCPHO',wscpho,0.0d0)
+      call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
 
 ! 12/1/95 Added weight for the multi-body term WCORR
       call reada(weightcard,'WCORRH',wcorr,1.0D0)

diff --git a/source/wham/conform_compar.f90 b/source/wham/conform_compar.f90
index 0470c67..8476646 100644 (file)
--- a/source/wham/conform_compar.f90
+++ b/source/wham/conform_compar.f90
@@ -1012,7 +1012,7 @@
 !----------------------------------------------------------------------------
       subroutine contact(lprint,ncont,icont)
 
-      use energy_data, only: nnt,nct,itype,ipot,maxcont,sigma,sigmaii
+      use energy_data, only: nnt,nct,itype,ipot,maxcont,sigma,sigmaii,molnum
 !      include 'DIMENSIONS'
 !      include 'COMMON.IOUNITS'
 !      include 'COMMON.CHAIN'
@@ -1024,6 +1024,7 @@
       logical :: lprint
       integer :: kkk,i,j,i1,i2,it1,it2,iti,itj
       real(kind=8) :: rcomp
+      integer :: mnum,mnum2
       ncont=0
       kkk=3
 !     print *,'nnt=',nnt,' nct=',nct
@@ -3178,8 +3179,8 @@
 !      include 'COMMON.CHAIN'
 !      include 'COMMON.NAMES'
 !      include 'COMMON.INTERACT'
-      integer :: ncont,icont(2,maxcont),isec(nres,4),nsec(nres),&
-        isecstr(nres)
+      integer :: ncont,icont(2,maxcont),isec(nres,4),nsec(nres+1),&
+        isecstr(nres+1)
       logical :: lprint,lprint_sec,not_done !el,freeres
       integer :: i,j,ii1,jj1,i1,j1,ij,k,ien,ist
       integer :: nstrand,nbeta,nhelix,iii1,jjj1,mnum

diff --git a/source/wham/enecalc.f90 b/source/wham/enecalc.f90
index 597b7cc..b8fd48f 100644 (file)
--- a/source/wham/enecalc.f90
+++ b/source/wham/enecalc.f90
@@ -3,7 +3,7 @@
       use io_units
       use wham_data
 !
-      use geometry_data, only:nres
+      use geometry_data, only:nres,boxxsize,boxysize,boxzsize
       use energy_data
 !      use control_data, only:maxthetyp1
       use energy, only:etotal,enerprint,rescale_weights
@@ -95,7 +95,7 @@
       use geometry_data, only:c,phi,theta,alph,omeg,deg2rad,anatemp,&
                               vbld,rad2deg,dc_norm,dc,vbld_inv
       use io_base, only:gyrate!,briefout
-      use geometry, only:int_from_cart1
+      use geometry, only:int_from_cart1,returnbox
       use io_wham, only:pdboutW
       use io_database, only:opentmp
       use conform_compar, only:qwolynes,rmsnat
@@ -211,6 +211,17 @@ write(iout,*)"enecalc_ i ntot",i,ntot
                c(l,k+nres)=csingle(l,k+nres)
              enddo
            enddo
+           call returnbox
+           do k=1,nres
+             do l=1,3
+               csingle(l,k)=c(l,k)
+             enddo
+           enddo
+           do k=nnt,nct
+             do l=1,3
+               csingle(l,k+nres)=c(l,k+nres)
+             enddo
+           enddo
            anatemp= 1.0d0/(beta_h(ib,ipar)*1.987D-3)
            q(nQ+1,iii+1)=rmsnat(iii+1)
          endif
@@ -1169,6 +1180,7 @@ write(iout,*)"end of store_parm"
       use geometry_data, only:c
       use io_base, only:ilen
       use io_wham, only:pdboutW
+      use geometry, only:returnbox
 !      implicit none
 !      include "DIMENSIONS"
 !      include "DIMENSIONS.ZSCOPT"
@@ -1344,6 +1356,9 @@ write(iout,*)"end of store_parm"
             esccor=enetb(21,i,iparm)
 !            edihcnstr=enetb(20,i,iparm)
             edihcnstr=enetb(19,i,iparm)
+            ecationcation=enetb(42,i,iparm)
+            ecation_prot=enetb(41,i,iparm)
+
 !#ifdef SPLITELE
 !            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
 !            +wvdwpp*evdw1 &
@@ -1542,6 +1557,13 @@ write(iout,*)"end of store_parm"
                 c(k,j)=csingle(k,j)
               enddo
             enddo
+            call returnbox
+            do j=1,2*nres
+              do k=1,3
+                csingle(k,j)=c(k,j)
+              enddo
+            enddo
+ 
             eini=fdimless(i)
             call pdboutW(iperm(i),temper,eini,enepot(i),efree,rmsdev)
           enddo