unres_package_Oct_2016 from emilial
[unres4.git] / source / wham / enecalc.f90
diff --git a/source/wham/enecalc.f90 b/source/wham/enecalc.f90
new file mode 100644 (file)
index 0000000..fd5f6ca
--- /dev/null
@@ -0,0 +1,1708 @@
+      module ene_calc
+!-----------------------------------------------------------------------------
+      use io_units
+      use wham_data
+!
+      use geometry_data, only:nres
+      use energy_data
+      use control_data, only:maxthetyp1
+      use energy, only:etotal,enerprint,rescale_weights
+#ifdef MPI
+      use MPI_data
+!      include "mpif.h"
+!      include "COMMON.MPI"
+#endif
+      implicit none
+!-----------------------------------------------------------------------------
+! COMMON.ALLPARM
+!      common /allparm/
+      real(kind=8),dimension(:,:),allocatable :: ww_all !(max_ene,max_parm) ! max_eneW
+      real(kind=8),dimension(:),allocatable :: vbldp0_all,akp_all !(max_parm)
+      real(kind=8),dimension(:,:,:),allocatable :: vbldsc0_all,&
+        aksc_all,abond0_all !(maxbondterm,ntyp,max_parm)
+      real(kind=8),dimension(:,:),allocatable :: a0thet_all !(-ntyp:ntyp,max_parm)
+      real(kind=8),dimension(:,:,:,:,:),allocatable :: athet_all,&
+        bthet_all !(2,-ntyp:ntyp,-1:1,-1:1,max_parm)
+      real(kind=8),dimension(:,:,:),allocatable :: polthet_all !(0:3,-ntyp:ntyp,max_parm)
+      real(kind=8),dimension(:,:,:),allocatable :: gthet_all !(3,-ntyp:ntyp,max_parm)
+      real(kind=8),dimension(:,:),allocatable :: theta0_all,&
+        sig0_all,sigc0_all !(-ntyp:ntyp,max_parm)
+      real(kind=8),dimension(:,:,:,:,:),allocatable :: aa0thet_all
+!(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
+      real(kind=8),dimension(:,:,:,:,:,:),allocatable :: aathet_all
+!(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
+      real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: bbthet_all,&
+        ccthet_all,ddthet_all,eethet_all !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,
+!     & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
+      real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: ffthet_all1,&
+        ggthet_all1,ffthet_all2,ggthet_all2 !(maxdouble,maxdouble,maxtheterm3,
+!     &  -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm)
+      real(kind=8),dimension(:,:),allocatable :: dsc_all,dsc0_all !(ntyp1,max_parm)
+      real(kind=8),dimension(:,:,:),allocatable :: bsc_all !(maxlob,ntyp,max_parm)
+      real(kind=8),dimension(:,:,:,:),allocatable :: censc_all !(3,maxlob,-ntyp:ntyp,max_parm)
+      real(kind=8),dimension(:,:,:,:,:),allocatable :: gaussc_all !(3,3,maxlob,-ntyp:ntyp,max_parm)
+      real(kind=8),dimension(:,:,:),allocatable :: sc_parmin_all !(65,ntyp,max_parm)
+      real(kind=8),dimension(:,:,:,:),allocatable :: v0_all
+!(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
+      real(kind=8),dimension(:,:,:,:,:),allocatable :: v1_all,&
+        v2_all !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
+      real(kind=8),dimension(:,:,:,:),allocatable :: vlor1_all,&
+        vlor2_all,vlor3_all !(maxlor,maxtor,maxtor,max_parm)
+      real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v1c_all,&
+        v1s_all !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
+      real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v2c_all
+!(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
+      real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v2s_all
+!(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
+      real(kind=8),dimension(:,:,:),allocatable :: b1_all,b2_all !(2,-maxtor:maxtor,max_parm)
+      real(kind=8),dimension(:,:,:,:),allocatable :: cc_all,dd_all,&
+        ee_all !(2,2,-maxtor:maxtor,max_parm)
+      real(kind=8),dimension(:,:,:,:),allocatable :: ctilde_all,&
+        dtilde_all !(2,2,-maxtor:maxtor,max_parm)
+      real(kind=8),dimension(:,:,:),allocatable :: b1tilde_all !(2,-maxtor:maxtor,max_parm)
+      real(kind=8),dimension(:,:,:),allocatable :: app_all,bpp_all,&
+        ael6_all,ael3_all !(2,2,max_parm)
+      real(kind=8),dimension(:,:,:),allocatable :: aad_all,&
+        bad_all !(ntyp,2,max_parm)
+      real(kind=8),dimension(:,:,:),allocatable :: aa_all,bb_all,&
+        augm_all,eps_all,sigma_all,r0_all,chi_all !(ntyp,ntyp,max_parm)
+      real(kind=8),dimension(:,:),allocatable :: chip_all,alp_all !(ntyp,max_parm)
+      real(kind=8),dimension(:),allocatable :: ebr_all,d0cm_all,&
+        akcm_all,akth_all,akct_all,v1ss_all,v2ss_all,v3ss_all !(max_parm)
+      real(kind=8),dimension(:,:,:,:,:),allocatable :: v1sccor_all,&
+        v2sccor_all !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
+      integer,dimension(:,:),allocatable :: nlob_all !(ntyp1,max_parm)
+      integer,dimension(:,:,:,:),allocatable :: nlor_all,&
+        nterm_all !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
+      integer,dimension(:,:,:,:,:),allocatable :: ntermd1_all,&
+        ntermd2_all !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
+      integer,dimension(:,:),allocatable :: nbondterm_all !(ntyp,max_parm)
+      integer,dimension(:,:),allocatable :: ithetyp_all !(-ntyp1:ntyp1,max_parm)
+      integer,dimension(:),allocatable :: nthetyp_all,ntheterm_all,&
+        ntheterm2_all,ntheterm3_all,nsingle_all,ndouble_all,&
+        nntheterm_all !(max_parm)
+      integer,dimension(:,:,:),allocatable :: nterm_sccor_all !(-ntyp:ntyp,-ntyp:ntyp,max_parm)
+!-----------------------------------------------------------------------------
+!
+!
+!-----------------------------------------------------------------------------
+      contains
+!-----------------------------------------------------------------------------
+      subroutine enecalc(islice,*)
+
+      use names
+      use control_data, only:indpdb
+      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 io_wham, only:pdboutW
+      use io_database, only:opentmp
+      use conform_compar, only:qwolynes,rmsnat
+!      include "DIMENSIONS"
+!      include "DIMENSIONS.ZSCOPT"
+!      include "DIMENSIONS.FREE"
+#ifdef MPI
+!      use MPI_data
+      include "mpif.h"
+!      include "COMMON.MPI"
+#endif
+!      include "COMMON.CHAIN"
+!      include "COMMON.IOUNITS"
+!      include "COMMON.PROTFILES"
+!      include "COMMON.NAMES"
+!      include "COMMON.VAR"
+!      include "COMMON.SBRIDGE"
+!      include "COMMON.GEO"
+!      include "COMMON.FFIELD"
+!      include "COMMON.ENEPS"
+!      include "COMMON.LOCAL"
+!      include "COMMON.WEIGHTS"
+!      include "COMMON.INTERACT"
+!      include "COMMON.FREE"
+!      include "COMMON.ENERGIES"
+!      include "COMMON.CONTROL"
+!      include "COMMON.TORCNSTR"
+!      implicit none
+#ifdef MPI
+      integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
+#endif
+      character(len=64) :: nazwa
+      character(len=80) :: bxname
+      character(len=3) :: liczba
+!el      real(kind=8) :: qwolynes
+!el      external qwolynes
+      integer :: errmsg_count,maxerrmsg_count=100
+!el      real(kind=8) :: rmsnat,gyrate
+!el      external rmsnat,gyrate
+      real(kind=8) :: tole=1.0d-1
+      integer i,itj,ii,iii,j,k,l,licz
+      integer ir,ib,ipar,iparm
+      integer iscor,islice
+      real(kind=4) :: csingle(3,nres*2)
+      real(kind=8) :: energ
+      real(kind=8) :: temp
+!el      integer ilen,iroof
+!el      external ilen,iroof
+      real(kind=8) :: energia(0:n_ene),rmsdev,efree,eini
+!el      real(kind=8) :: energia(0:max_ene),rmsdev,efree,eini
+      real(kind=8) :: fT(6),quot,quotl,kfacl,kfac=2.4d0,T0=3.0d2
+      real(kind=8) :: tt
+      integer :: snk_p(MaxR,MaxT_h,nParmSet)!Max_parm)
+      logical :: lerr
+      character(len=64) :: bprotfile_temp
+
+!      integer :: rec
+      integer,dimension(0:nprocs) :: scount_
+!el      real(kind=8) :: rmsnat
+
+      rescale_mode=rescale_modeW
+
+      call opentmp(islice,ientout,bprotfile_temp)
+      iii=0
+      ii=0
+!el
+!      iparm=1
+      errmsg_count=0
+      write (iout,*) "enecalc: nparmset ",nparmset
+#ifdef MPI
+      do iparm=1,nParmSet
+        do ib=1,nT_h(iparm)
+          do i=1,nR(ib,iparm)
+            snk_p(i,ib,iparm)=0
+          enddo
+        enddo
+      enddo
+      do i=indstart(me1),indend(me1)
+write(iout,*)"enecalc_ i indstart",i,indstart(me1),indend(me1)
+#else
+      do iparm=1,nParmSet
+        do ib=1,nT_h(iparm)
+          do i=1,nR(ib,iparm)
+            snk(i,ib,iparm)=0
+          enddo
+        enddo
+      enddo
+      do i=1,ntot
+write(iout,*)"enecalc_ i ntot",i,ntot
+#endif
+        read(ientout,rec=i,err=101) &
+          ((csingle(l,k),l=1,3),k=1,nres),&
+          ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
+          nss,(ihpb(k),jhpb(k),k=1,nss),&
+          eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
+!el debug
+!write(iout,*)"co wczytuje"
+!          write(iout,*)((csingle(l,k),l=1,3),k=1,nres),&
+!          ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
+!          nss,(ihpb(k),jhpb(k),k=1,nss),&
+!          eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
+!el --------
+
+!write(iout,*)"ipar",ib,ipar,1.0d0/(beta_h(ib,ipar)*1.987D-3)
+         if (indpdb.gt.0) then
+           do k=1,nres
+             do l=1,3
+               c(l,k)=csingle(l,k)
+             enddo
+           enddo
+           do k=nnt,nct
+             do l=1,3
+               c(l,k+nres)=csingle(l,k+nres)
+             enddo
+           enddo
+           anatemp= 1.0d0/(beta_h(ib,ipar)*1.987D-3)
+           q(nQ+1,iii+1)=rmsnat(iii+1)
+         endif
+         q(nQ+2,iii+1)=gyrate(iii+1)
+! write(iout,*)"wczyt",anatemp,q(nQ+2,iii+1) !el
+!        fT=T0*beta_h(ib,ipar)*1.987D-3
+!        ft=2.0d0/(1.0d0+1.0d0/(T0*beta_h(ib,ipar)*1.987D-3))
+! EL start old rescale
+!        if (rescale_modeW.eq.1) then
+!          quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
+!#if defined(FUNCTH)
+!          tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
+!          ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
+!#elif defined(FUNCT)
+!          ft(6)=quot
+!#else
+!          ft(6)=1.0d0
+!#endif
+!          quotl=1.0d0
+!          kfacl=1.0d0
+!          do l=1,5
+!            quotl=quotl*quot
+!            kfacl=kfacl*kfac
+!            fT(l)=kfacl/(kfacl-1.0d0+quotl)
+!          enddo
+!        else if (rescale_modeW.eq.2) then
+!          quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
+!#if defined(FUNCTH)
+!          tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
+!          ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
+!#elif defined(FUNCT)
+!          ft(6)=quot
+!#else
+!          ft(6)=1.0d0
+!#endif
+!          quotl=1.0d0
+!          do l=1,5
+!            quotl=quotl*quot
+!            fT(l)=1.12692801104297249644d0/ &
+!               dlog(dexp(quotl)+dexp(-quotl))
+!          enddo
+!        else if (rescale_modeW.eq.0) then
+!          do l=1,5
+!            fT(l)=1.0d0
+!          enddo
+!        else
+!          write (iout,*) "Error in ECECALC: wrong RESCALE_MODE",&
+!           rescale_modeW
+!          call flush(iout)
+!          return 1
+!        endif
+!EL end old rescele
+!        write (iout,*) "T",1.0d0/(beta_h(ib,ipar)*1.987D-3)," T0",T0,
+!     &   " kfac",kfac,"quot",quot," fT",fT
+#ifdef DEBUG
+            write(iout,*)"weights"
+            write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
+            wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
+            wtor_d,wsccor,wbond
+#endif
+
+        do j=1,2*nres
+          do k=1,3
+            c(k,j)=csingle(k,j)
+          enddo
+        enddo
+        call int_from_cart1(.false.)
+        ii=ii+1
+
+!        call rescale_weights(1.0d0/(beta_h(ib,ipar)*1.987D-3))
+       do iparm=1,nparmset
+#ifdef DEBUG
+            write (iout,*) "before restore w=",1.0d0/(beta_h(ib,ipar)*1.987D-3)
+            write(iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
+            wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
+            wtor_d,wsccor,wbond
+#endif
+        call restore_parm(iparm)
+        call rescale_weights(1.0d0/(beta_h(ib,ipar)*1.987D-3))
+#ifdef DEBUG
+            write (iout,*) "before etot w=",1.0d0/(beta_h(ib,ipar)*1.987D-3)
+            write(iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
+            wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
+            wtor_d,wsccor,wbond
+#endif
+!        call etotal(energia(0),fT)
+        call etotal(energia(0))
+!write(iout,*)"check c and dc after etotal",1.0d0/(0.001987*beta_h(ib,ipar))
+!do k=1,2*nres+2
+!write(iout,*)k,"c=",(c(l,k),l=1,3)
+!write(iout,*)k,"dc=",(dc(l,k),l=1,3)
+!write(iout,*)k,"dc_norm=",(dc_norm(l,k),l=1,3)
+!enddo
+!do k=1,nres*2
+!write(iout,*)k,"vbld=",vbld(k)
+!write(iout,*)k,"vbld_inv=",vbld_inv(k)
+!enddo
+
+!write(iout,*)"energia",(energia(j),j=0,n_ene)
+!write(iout,*)"enerprint tuz po call etotal"
+        call enerprint(energia(0))
+#ifdef DEBUG
+        write (iout,*) "Conformation",i
+          write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
+          write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
+!        call enerprint(energia(0),fT)
+        call enerprint(energia(0))
+        write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,21)
+        write (iout,*) "ftors",ftors
+!el        call briefout(i,energia(0))
+        temp=1.0d0/(beta_h(ib,ipar)*1.987D-3)
+        write (iout,*) "temp", temp
+        call pdboutW(i,temp,energia(0),energia(0),0.0d0,0.0d0)
+#endif
+        if (energia(0).ge.1.0d20) then
+          write (iout,*) "NaNs detected in some of the energy",&
+           " components for conformation",ii+1
+          write (iout,*) "The Cartesian geometry is:"
+          write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
+          write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
+          write (iout,*) "The internal geometry is:"
+!          call intout
+!        call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
+          write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
+          write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
+          write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
+          write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
+          write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
+          write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
+          write (iout,*) "The components of the energy are:"
+!          call enerprint(energia(0),fT)
+          call enerprint(energia(0))
+          write (iout,*) &
+            "This conformation WILL NOT be added to the database."
+          call flush(iout)
+          goto 121
+        else 
+#ifdef DEBUG
+          if (ipar.eq.iparm) write (iout,*) i,iparm,&
+            1.0d0/(beta_h(ib,ipar)*1.987D-3),eini,energia(0)
+#endif
+          if (ipar.eq.iparm .and. einicheck.gt.0 .and. &
+            dabs(eini-energia(0)).gt.tole) then
+            if (errmsg_count.le.maxerrmsg_count) then
+              write (iout,'(2a,2e15.5,a,2i8,a,f8.1)') &
+               "Warning: energy differs remarkably from ",&
+               " the value read in: ",energia(0),eini," point",&
+               iii+1,indstart(me1)+iii," T",&
+               1.0d0/(1.987D-3*beta_h(ib,ipar))
+!              call intout
+              call pdboutW(indstart(me1)+iii,&
+       1.0d0/(1.987D-3*beta_h(ib,ipar)),&
+       energia(0),eini,0.0d0,0.0d0)
+          write (iout,*) "The Cartesian geometry is:"
+          write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
+          write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
+          write (iout,*) "The internal geometry is:"
+!          call intout
+!        call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
+          write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
+          write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
+          write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
+          write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
+          write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
+          write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
+              call enerprint(energia(0))
+!              call enerprint(energia(0),fT)
+              errmsg_count=errmsg_count+1
+              if (errmsg_count.gt.maxerrmsg_count) &
+                write (iout,*) "Too many warning messages"
+              if (einicheck.gt.1) then
+                write (iout,*) "Calculation stopped."
+                call flush(iout)
+#ifdef MPI
+                call MPI_Abort(WHAM_COMM,IERROR,ERRCODE)
+#endif
+                call flush(iout)
+                return 1
+              endif
+            endif
+          endif
+          potE(iii+1,iparm)=energia(0)
+          do k=1,21
+            enetb(k,iii+1,iparm)=energia(k)
+          enddo
+!           write (iout,'(2i5,21f8.2)') "debug",k,iii+1,(enetb(k,iii+1,iparm),k=1,21)
+!           write (iout,*) "debug",k,iii+1,(enetb(k,iii+1,iparm),k=1,21)
+#ifdef DEBUG
+          write (iout,'(2i5,f10.1,3e15.5)') i,iii,&
+           1.0d0/(beta_h(ib,ipar)*1.987D-3),energia(0),eini,efree
+!          call enerprint(energia(0),fT)
+#endif
+#ifdef DEBUG
+          write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
+          write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
+          write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
+          write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
+          write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
+          write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
+          write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
+          write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
+          write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
+          write (iout,'(8f10.5)') (q(k,iii+1),k=1,nQ)
+          write (iout,'(f10.5,i10)') rmsdev,iscor
+!          call enerprint(energia(0),fT)
+          call enerprint(energia(0))
+        call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
+#endif
+        endif
+
+        enddo ! iparm
+
+        iii=iii+1
+        if (q(1,iii).le.0.0d0 .and. indpdb.gt.0) q(1,iii)=qwolynes(0,0)
+        write (ientout,rec=iii) &
+         ((csingle(l,k),l=1,3),k=1,nres),&
+         ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
+         nss,(ihpb(k),jhpb(k),k=1,nss),&
+         potE(iii,ipar),efree,rmsdev,(q(k,iii),k=1,nQ),iR,ib,ipar
+!        write (iout,'(2i5,2e15.5)') ii,iii,potE(iii,ipar),efree
+#ifdef MPI
+        if (separate_parset) then
+          snk_p(iR,ib,1)=snk_p(iR,ib,1)+1
+        else
+          snk_p(iR,ib,ipar)=snk_p(iR,ib,ipar)+1
+        endif
+!        write (iout,*) "iii",iii," iR",iR," ib",ib," ipar",ipar,
+!     &   " snk",snk_p(iR,ib,ipar)
+#else
+        snk(iR,ib,ipar,islice)=snk(iR,ib,ipar,islice)+1
+#endif
+  121   continue
+      enddo   
+#ifdef MPI
+      scount(me)=iii 
+      write (iout,*) "Me",me," scount",scount(me)
+      call flush(iout)
+!  Master gathers updated numbers of conformations written by all procs.
+      call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount_(0), 1, &
+        MPI_INTEGER, WHAM_COMM, IERROR)
+      indstart(0)=1
+      indend(0)=scount_(0)
+      do i=1, Nprocs-1
+        indstart(i)=indend(i-1)+1
+        indend(i)=indstart(i)+scount_(i)-1
+      enddo
+      write (iout,*)
+      write (iout,*) "Revised conformation counts"
+      do i=0,nprocs1-1
+        write (iout,'(a,i5,a,i7,a,i7,a,i7)') &
+          "Processor",i," indstart",indstart(i),&
+          " indend",indend(i)," count",scount_(i)
+      enddo
+      call flush(iout)
+      call MPI_AllReduce(snk_p(1,1,1),snk(1,1,1,islice),&
+        MaxR*MaxT_h*nParmSet,&
+        MPI_INTEGER,MPI_SUM,WHAM_COMM,IERROR)
+#endif
+      stot(islice)=0
+      do iparm=1,nParmSet
+        do ib=1,nT_h(iparm)
+          do i=1,nR(ib,iparm)
+            stot(islice)=stot(islice)+snk(i,ib,iparm,islice)
+          enddo
+        enddo
+      enddo
+      write (iout,*) "Revised SNK"
+      do iparm=1,nParmSet
+        do ib=1,nT_h(iparm)
+          write (iout,'("Param",i3," Temp",f6.1,3x,32i8)') &
+           iparm,1.0d0/(1.987D-3*beta_h(ib,iparm)),&
+           (snk(i,ib,iparm,islice),i=1,nR(ib,iparm))
+          write (iout,*) "snk_p",(snk_p(i,ib,iparm),i=1,nR(ib,iparm))
+        enddo
+      enddo
+      write (iout,'("Total",i10)') stot(islice)
+      call flush(iout)
+      do i=0,nprocs
+        scount(i)=scount_(i)
+      enddo
+      return
+  101 write (iout,*) "Error in scratchfile."
+      call flush(iout)
+!el#undef DEBUG
+      return 1
+      end subroutine enecalc
+!------------------------------------------------------------------------------
+      logical function conf_check(ii,iprint)
+
+      use geometry_data, only:c,phi,theta,alph,omeg,deg2rad,rad2deg,vbld
+      use geometry, only:int_from_cart1
+!      include "DIMENSIONS"
+!      include "DIMENSIONS.ZSCOPT"
+!      include "DIMENSIONS.FREE"
+#ifdef MPI
+!      use MPI_data
+      include "mpif.h"
+!      include "COMMON.MPI"
+#endif
+!      include "COMMON.CHAIN"
+!      include "COMMON.IOUNITS"
+!      include "COMMON.PROTFILES"
+!      include "COMMON.NAMES"
+!      include "COMMON.VAR"
+!      include "COMMON.SBRIDGE"
+!      include "COMMON.GEO"
+!      include "COMMON.FFIELD"
+!      include "COMMON.ENEPS"
+!      include "COMMON.LOCAL"
+!      include "COMMON.WEIGHTS"
+!      include "COMMON.INTERACT"
+!      include "COMMON.FREE"
+!      include "COMMON.ENERGIES"
+!      include "COMMON.CONTROL"
+!      include "COMMON.TORCNSTR"
+!      implicit none
+#ifdef MPI
+      integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
+#endif
+      integer :: j,k,l,ii,itj,iprint
+      if (.not.check_conf) then
+        conf_check=.true.
+        return
+      endif
+      call int_from_cart1(.false.)
+      do j=nnt+1,nct
+        if (itype(j-1).ne.ntyp1 .and. itype(j).ne.ntyp1 .and. &
+          (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
+          if (iprint.gt.0) &
+          write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),&
+            " for conformation",ii
+          if (iprint.gt.1) then
+            write (iout,*) "The Cartesian geometry is:"
+            write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
+            write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
+            write (iout,*) "The internal geometry is:"
+            write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
+            write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
+            write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
+            write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
+            write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
+            write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
+          endif
+          if (iprint.gt.0) write (iout,*) &
+            "This conformation WILL NOT be added to the database."
+          conf_check=.false.
+          return
+        endif
+      enddo
+      do j=nnt,nct
+        itj=itype(j)
+        if (itype(j).ne.10 .and.itype(j).ne.ntyp1 .and. &
+           (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then
+          if (iprint.gt.0) &
+          write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),&
+           " for conformation",ii
+          if (iprint.gt.1) then
+            write (iout,*) "The Cartesian geometry is:"
+            write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
+            write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
+            write (iout,*) "The internal geometry is:"
+            write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
+            write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
+            write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
+            write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
+            write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
+            write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
+          endif
+          if (iprint.gt.0) write (iout,*) &
+            "This conformation WILL NOT be added to the database."
+          conf_check=.false.
+          return
+        endif
+      enddo
+      do j=3,nres
+        if (theta(j).le.0.0d0) then
+          if (iprint.gt.0) &
+          write (iout,*) "Zero theta angle(s) in conformation",ii
+          if (iprint.gt.1) then
+            write (iout,*) "The Cartesian geometry is:"
+            write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
+            write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
+            write (iout,*) "The internal geometry is:"
+            write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
+            write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
+            write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
+            write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
+            write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
+            write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
+          endif
+          if (iprint.gt.0) write (iout,*) &
+            "This conformation WILL NOT be added to the database." 
+          conf_check=.false.
+          return
+        endif
+        if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
+      enddo
+      conf_check=.true.
+!      write (iout,*) "conf_check passed",ii
+      return
+      end function conf_check
+!-----------------------------------------------------------------------------
+! store_parm.F
+!-----------------------------------------------------------------------------
+      subroutine store_parm(iparm)
+!
+! Store parameters of set IPARM
+! valence angles and the side chains and energy parameters.
+!
+!      implicit none
+!      include 'DIMENSIONS'
+!      include 'DIMENSIONS.ZSCOPT'
+!      include 'DIMENSIONS.FREE'
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.CHAIN'
+!      include 'COMMON.INTERACT'
+!      include 'COMMON.GEO'
+!      include 'COMMON.LOCAL'
+!      include 'COMMON.TORSION'
+!      include 'COMMON.FFIELD'
+!      include 'COMMON.NAMES'
+!      include 'COMMON.SBRIDGE'
+!      include 'COMMON.SCROT'
+!      include 'COMMON.SCCOR'
+!      include 'COMMON.ALLPARM'
+      integer :: i,j,k,l,m,mm,iparm,ichir1,ichir2,iblock,iii
+
+      call alloc_enecalc_arrays(iparm)
+!el      allocate(ww_all(n_ene,iparm))
+! Store weights
+      ww_all(1,iparm)=wsc
+      ww_all(2,iparm)=wscp
+      ww_all(3,iparm)=welec
+      ww_all(4,iparm)=wcorr
+      ww_all(5,iparm)=wcorr5
+      ww_all(6,iparm)=wcorr6
+      ww_all(7,iparm)=wel_loc
+      ww_all(8,iparm)=wturn3
+      ww_all(9,iparm)=wturn4
+      ww_all(10,iparm)=wturn6
+      ww_all(11,iparm)=wang
+      ww_all(12,iparm)=wscloc
+      ww_all(13,iparm)=wtor
+      ww_all(14,iparm)=wtor_d
+      ww_all(15,iparm)=wstrain
+      ww_all(16,iparm)=wvdwpp
+      ww_all(17,iparm)=wbond
+      ww_all(19,iparm)=wsccor
+! Store bond parameters
+      vbldp0_all(iparm)=vbldp0
+      akp_all(iparm)=akp
+      do i=1,ntyp
+        nbondterm_all(i,iparm)=nbondterm(i)
+        do j=1,nbondterm(i)
+          vbldsc0_all(j,i,iparm)=vbldsc0(j,i)
+          aksc_all(j,i,iparm)=aksc(j,i)
+          abond0_all(j,i,iparm)=abond0(j,i)
+        enddo
+      enddo
+! Store bond angle parameters
+#ifdef CRYST_THETA
+      do i=-ntyp,ntyp
+        a0thet_all(i,iparm)=a0thet(i)
+        do ichir1=-1,1
+        do ichir2=-1,1
+        do j=1,2
+          athet_all(j,i,ichir1,ichir2,iparm)=athet(j,i,ichir1,ichir2)
+          bthet_all(j,i,ichir1,ichir2,iparm)=bthet(j,i,ichir1,ichir2)
+        enddo
+        enddo
+        enddo
+        do j=0,3
+          polthet_all(j,i,iparm)=polthet(j,i)
+        enddo
+        do j=1,3
+          gthet_all(j,i,iparm)=gthet(j,i)
+        enddo
+        theta0_all(i,iparm)=theta0(i)
+        sig0_all(i,iparm)=sig0(i)
+        sigc0_all(i,iparm)=sigc0(i)
+      enddo
+#else
+      nthetyp_all(iparm)=nthetyp
+      ntheterm_all(iparm)=ntheterm
+      ntheterm2_all(iparm)=ntheterm2
+      ntheterm3_all(iparm)=ntheterm3
+      nsingle_all(iparm)=nsingle
+      ndouble_all(iparm)=ndouble
+      nntheterm_all(iparm)=nntheterm
+      do i=-ntyp,ntyp
+        ithetyp_all(i,iparm)=ithetyp(i)
+      enddo
+      do iblock=1,2
+      do i=-maxthetyp1,maxthetyp1
+        do j=-maxthetyp1,maxthetyp1
+          do k=-maxthetyp1,maxthetyp1
+            aa0thet_all(i,j,k,iblock,iparm)=aa0thet(i,j,k,iblock)
+            do l=1,ntheterm
+              aathet_all(l,i,j,k,iblock,iparm)=aathet(l,i,j,k,iblock)
+            enddo
+            do l=1,ntheterm2
+              do m=1,nsingle
+                bbthet_all(m,l,i,j,k,iblock,iparm)= &
+      bbthet(m,l,i,j,k,iblock)
+                ccthet_all(m,l,i,j,k,iblock,iparm)= &
+      ccthet(m,l,i,j,k,iblock)
+                ddthet_all(m,l,i,j,k,iblock,iparm)= &
+      ddthet(m,l,i,j,k,iblock)
+                eethet_all(m,l,i,j,k,iblock,iparm)= &
+      eethet(m,l,i,j,k,iblock)
+              enddo
+            enddo
+            do l=1,ntheterm3
+              do m=1,ndouble
+                do mm=1,ndouble
+                if (iblock.eq.1) then
+                 ffthet_all1(mm,m,l,i,j,k,iparm)=&
+         ffthet(mm,m,l,i,j,k,iblock)
+                 ggthet_all1(mm,m,l,i,j,k,iparm)=&
+      ggthet(mm,m,l,i,j,k,iblock)
+                  else
+                 ffthet_all2(mm,m,l,i,j,k,iparm)=&
+         ffthet(mm,m,l,i,j,k,iblock)
+                 ggthet_all2(mm,m,l,i,j,k,iparm)=&
+      ggthet(mm,m,l,i,j,k,iblock)
+                  endif
+                enddo
+              enddo
+            enddo
+          enddo
+        enddo
+      enddo
+      enddo
+#endif
+#ifdef CRYST_SC
+! Store the sidechain rotamer parameters
+      do i=-ntyp,ntyp
+       iii=iabs(i)
+!!       write (iout,*) i,"storeparm1"
+       if (i.eq.0) cycle
+        nlob_all(iii,iparm)=nlob(iii)
+        do j=1,nlob(iii)
+          bsc_all(j,iii,iparm)=bsc(j,iii)
+          do k=1,3
+            censc_all(k,j,i,iparm)=censc(k,j,i)
+          enddo
+          do k=1,3
+            do l=1,3
+              gaussc_all(l,k,j,i,iparm)=gaussc(l,k,j,i)
+            enddo
+          enddo
+        enddo
+      enddo
+#else
+      do i=1,ntyp
+        do j=1,65
+          sc_parmin_all(j,i,iparm)=sc_parmin(j,i)
+        enddo
+      enddo
+#endif
+! Store the torsional parameters
+      do iblock=1,2
+      do i=-ntortyp+1,ntortyp-1
+        do j=-ntortyp+1,ntortyp-1
+          v0_all(i,j,iblock,iparm)=v0(i,j,iblock)
+          nterm_all(i,j,iblock,iparm)=nterm(i,j,iblock)
+          nlor_all(i,j,iblock,iparm)=nlor(i,j,iblock)
+          do k=1,nterm(i,j,iblock)
+            v1_all(k,i,j,iblock,iparm)=v1(k,i,j,iblock)
+            v2_all(k,i,j,iblock,iparm)=v2(k,i,j,iblock)
+          enddo
+          do k=1,nlor(i,j,iblock)
+            vlor1_all(k,i,j,iparm)=vlor1(k,i,j)
+            vlor2_all(k,i,j,iparm)=vlor2(k,i,j)
+            vlor3_all(k,i,j,iparm)=vlor3(k,i,j)
+          enddo
+        enddo
+      enddo
+      enddo  
+! Store the double torsional parameters
+      do iblock=1,2
+      do i=-ntortyp+1,ntortyp-1
+        do j=-ntortyp+1,ntortyp-1
+          do k=-ntortyp+1,ntortyp-1
+            ntermd1_all(i,j,k,iblock,iparm)=ntermd_1(i,j,k,iblock)
+            ntermd2_all(i,j,k,iblock,iparm)=ntermd_2(i,j,k,iblock)
+            do l=1,ntermd_1(i,j,k,iblock)
+              v1c_all(1,l,i,j,k,iblock,iparm)=v1c(1,l,i,j,k,iblock)
+              v1c_all(2,l,i,j,k,iblock,iparm)=v1c(2,l,i,j,k,iblock)
+              v2c_all(1,l,i,j,k,iblock,iparm)=v2c(1,l,i,j,k,iblock)
+              v2c_all(2,l,i,j,k,iblock,iparm)=v2c(2,l,i,j,k,iblock)
+            enddo
+            do l=1,ntermd_2(i,j,k,iblock)
+              do m=1,ntermd_2(i,j,k,iblock)
+                v2s_all(l,m,i,j,k,iblock,iparm)=v2s(l,m,i,j,k,iblock)
+              enddo
+            enddo
+          enddo
+        enddo
+      enddo
+      enddo
+! Store parameters of the cumulants
+      do i=-nloctyp,nloctyp
+        do j=1,2
+          b1_all(j,i,iparm)=b1(j,i)
+          b1tilde_all(j,i,iparm)=b1tilde(j,i)
+          b2_all(j,i,iparm)=b2(j,i)
+        enddo
+        do j=1,2
+          do k=1,2
+            cc_all(k,j,i,iparm)=cc(k,j,i)
+            ctilde_all(k,j,i,iparm)=ctilde(k,j,i)
+            dd_all(k,j,i,iparm)=dd(k,j,i)
+            dtilde_all(k,j,i,iparm)=dtilde(k,j,i)
+            ee_all(k,j,i,iparm)=ee(k,j,i)
+          enddo
+        enddo
+      enddo
+! Store the parameters of electrostatic interactions
+      do i=1,2
+        do j=1,2
+          app_all(j,i,iparm)=app(j,i)
+          bpp_all(j,i,iparm)=bpp(j,i)
+          ael6_all(j,i,iparm)=ael6(j,i)
+          ael3_all(j,i,iparm)=ael3(j,i)
+        enddo
+      enddo
+! Store sidechain parameters
+      do i=1,ntyp
+        do j=1,ntyp
+          aa_all(j,i,iparm)=aa(j,i)
+          bb_all(j,i,iparm)=bb(j,i)
+          r0_all(j,i,iparm)=r0(j,i)
+          sigma_all(j,i,iparm)=sigma(j,i)
+          chi_all(j,i,iparm)=chi(j,i)
+          augm_all(j,i,iparm)=augm(j,i)
+          eps_all(j,i,iparm)=eps(j,i)
+        enddo
+      enddo
+      do i=1,ntyp
+        chip_all(i,iparm)=chip(i)
+        alp_all(i,iparm)=alp(i)
+      enddo
+! Store the SCp parameters
+      do i=1,ntyp
+        do j=1,2
+          aad_all(i,j,iparm)=aad(i,j)
+          bad_all(i,j,iparm)=bad(i,j)
+        enddo
+      enddo
+! Store disulfide-bond parameters
+      ebr_all(iparm)=ebr
+      d0cm_all(iparm)=d0cm
+      akcm_all(iparm)=akcm
+      akth_all(iparm)=akth
+      akct_all(iparm)=akct
+      v1ss_all(iparm)=v1ss
+      v2ss_all(iparm)=v2ss
+      v3ss_all(iparm)=v3ss
+! Store SC-backbone correlation parameters
+      do i=-nsccortyp,nsccortyp
+       do j=-nsccortyp,nsccortyp
+
+      nterm_sccor_all(j,i,iparm)=nterm_sccor(j,i)
+!      do i=1,20
+!        do j=1,20
+         do l=1,3
+          do k=1,nterm_sccor(j,i)
+            v1sccor_all(k,l,j,i,iparm)=v1sccor(k,l,j,i)
+            v2sccor_all(k,l,j,i,iparm)=v2sccor(k,l,j,i)
+           enddo
+          enddo
+        enddo
+      enddo
+write(iout,*)"end of store_parm"
+      return
+      end subroutine store_parm
+!--------------------------------------------------------------------------
+      subroutine restore_parm(iparm)
+!
+! Store parameters of set IPARM
+! valence angles and the side chains and energy parameters.
+!
+!      implicit none
+!      include 'DIMENSIONS'
+!      include 'DIMENSIONS.ZSCOPT'
+!      include 'DIMENSIONS.FREE'
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.CHAIN'
+!      include 'COMMON.INTERACT'
+!      include 'COMMON.GEO'
+!      include 'COMMON.LOCAL'
+!      include 'COMMON.TORSION'
+!      include 'COMMON.FFIELD'
+!      include 'COMMON.NAMES'
+!      include 'COMMON.SBRIDGE'
+!      include 'COMMON.SCROT'
+!      include 'COMMON.SCCOR'
+!      include 'COMMON.ALLPARM'
+      integer :: i,j,k,l,m,mm,iparm,ichir1,ichir2,iblock,iii
+
+! Restore weights
+      wsc=ww_all(1,iparm)
+      wscp=ww_all(2,iparm)
+      welec=ww_all(3,iparm)
+      wcorr=ww_all(4,iparm)
+      wcorr5=ww_all(5,iparm)
+      wcorr6=ww_all(6,iparm)
+      wel_loc=ww_all(7,iparm)
+      wturn3=ww_all(8,iparm)
+      wturn4=ww_all(9,iparm)
+      wturn6=ww_all(10,iparm)
+      wang=ww_all(11,iparm)
+      wscloc=ww_all(12,iparm)
+      wtor=ww_all(13,iparm)
+      wtor_d=ww_all(14,iparm)
+      wstrain=ww_all(15,iparm)
+      wvdwpp=ww_all(16,iparm)
+      wbond=ww_all(17,iparm)
+      wsccor=ww_all(19,iparm)
+! Restore bond parameters
+      vbldp0=vbldp0_all(iparm)
+      akp=akp_all(iparm)
+      do i=1,ntyp
+        nbondterm(i)=nbondterm_all(i,iparm)
+        do j=1,nbondterm(i)
+          vbldsc0(j,i)=vbldsc0_all(j,i,iparm)
+          aksc(j,i)=aksc_all(j,i,iparm)
+          abond0(j,i)=abond0_all(j,i,iparm)
+        enddo
+      enddo
+! Restore bond angle parameters
+#ifdef CRYST_THETA
+      do i=-ntyp,ntyp
+        a0thet(i)=a0thet_all(i,iparm)
+        do ichir1=-1,1
+        do ichir2=-1,1
+        do j=1,2
+          athet(j,i,ichir1,ichir2)=athet_all(j,i,ichir1,ichir2,iparm)
+          bthet(j,i,ichir1,ichir2)=bthet_all(j,i,ichir1,ichir2,iparm)
+        enddo
+        enddo
+        enddo
+        do j=0,3
+          polthet(j,i)=polthet_all(j,i,iparm)
+        enddo
+        do j=1,3
+          gthet(j,i)=gthet_all(j,i,iparm)
+        enddo
+        theta0(i)=theta0_all(i,iparm)
+        sig0(i)=sig0_all(i,iparm)
+        sigc0(i)=sigc0_all(i,iparm)
+      enddo
+#else
+      nthetyp=nthetyp_all(iparm)
+      ntheterm=ntheterm_all(iparm)
+      ntheterm2=ntheterm2_all(iparm)
+      ntheterm3=ntheterm3_all(iparm)
+      nsingle=nsingle_all(iparm)
+      ndouble=ndouble_all(iparm)
+      nntheterm=nntheterm_all(iparm)
+      do i=-ntyp,ntyp
+        ithetyp(i)=ithetyp_all(i,iparm)
+      enddo
+      do iblock=1,2
+      do i=-maxthetyp1,maxthetyp1
+        do j=-maxthetyp1,maxthetyp1
+          do k=-maxthetyp1,maxthetyp1
+            aa0thet(i,j,k,iblock)=aa0thet_all(i,j,k,iblock,iparm)
+            do l=1,ntheterm
+              aathet(l,i,j,k,iblock)=aathet_all(l,i,j,k,iblock,iparm)
+            enddo
+            do l=1,ntheterm2
+              do m=1,nsingle
+                bbthet(m,l,i,j,k,iblock)= &
+      bbthet_all(m,l,i,j,k,iblock,iparm)
+                ccthet(m,l,i,j,k,iblock)= &
+      ccthet_all(m,l,i,j,k,iblock,iparm)
+                ddthet(m,l,i,j,k,iblock)= &
+      ddthet_all(m,l,i,j,k,iblock,iparm)
+                eethet(m,l,i,j,k,iblock)= &
+      eethet_all(m,l,i,j,k,iblock,iparm)
+              enddo
+            enddo
+            do l=1,ntheterm3
+              do m=1,ndouble
+                do mm=1,ndouble
+                if (iblock.eq.1) then
+                 ffthet(mm,m,l,i,j,k,iblock)= &
+      ffthet_all1(mm,m,l,i,j,k,iparm)
+                 ggthet(mm,m,l,i,j,k,iblock)= &
+      ggthet_all1(mm,m,l,i,j,k,iparm)
+                else
+                 ffthet(mm,m,l,i,j,k,iblock)= &
+      ffthet_all2(mm,m,l,i,j,k,iparm)
+                 ggthet(mm,m,l,i,j,k,iblock)= &
+      ggthet_all2(mm,m,l,i,j,k,iparm)
+                endif
+                enddo
+              enddo
+            enddo
+          enddo
+        enddo
+      enddo
+      enddo
+#endif
+! Restore the sidechain rotamer parameters
+#ifdef CRYST_SC
+      do i=-ntyp,ntyp
+        if (i.eq.0) cycle
+        iii=iabs(i)
+        nlob(iii)=nlob_all(iii,iparm)
+        do j=1,nlob(iii)
+          bsc(j,iii)=bsc_all(j,iii,iparm)
+          do k=1,3
+            censc(k,j,i)=censc_all(k,j,i,iparm)
+          enddo
+          do k=1,3
+            do l=1,3
+              gaussc(l,k,j,i)=gaussc_all(l,k,j,i,iparm)
+            enddo
+          enddo
+        enddo
+      enddo
+#else
+      do i=1,ntyp
+        do j=1,65
+          sc_parmin(j,i)=sc_parmin_all(j,i,iparm)
+        enddo
+      enddo
+#endif
+! Restore the torsional parameters
+      do iblock=1,2
+      do i=-ntortyp+1,ntortyp-1
+        do j=-ntortyp+1,ntortyp-1
+          v0(i,j,iblock)=v0_all(i,j,iblock,iparm)
+          nterm(i,j,iblock)=nterm_all(i,j,iblock,iparm)
+          nlor(i,j,iblock)=nlor_all(i,j,iblock,iparm)
+          do k=1,nterm(i,j,iblock)
+            v1(k,i,j,iblock)=v1_all(k,i,j,iblock,iparm)
+            v2(k,i,j,iblock)=v2_all(k,i,j,iblock,iparm)
+          enddo
+          do k=1,nlor(i,j,iblock)
+            vlor1(k,i,j)=vlor1_all(k,i,j,iparm)
+            vlor2(k,i,j)=vlor2_all(k,i,j,iparm)
+            vlor3(k,i,j)=vlor3_all(k,i,j,iparm)
+          enddo
+        enddo
+      enddo  
+      enddo
+! Restore the double torsional parameters
+      do iblock=1,2
+      do i=-ntortyp+1,ntortyp-1
+        do j=-ntortyp+1,ntortyp-1
+          do k=-ntortyp+1,ntortyp-1
+            ntermd_1(i,j,k,iblock)=ntermd1_all(i,j,k,iblock,iparm)
+            ntermd_2(i,j,k,iblock)=ntermd2_all(i,j,k,iblock,iparm)
+            do l=1,ntermd_1(i,j,k,iblock)
+              v1c(1,l,i,j,k,iblock)=v1c_all(1,l,i,j,k,iblock,iparm)
+              v1c(2,l,i,j,k,iblock)=v1c_all(2,l,i,j,k,iblock,iparm)
+              v2c(1,l,i,j,k,iblock)=v2c_all(1,l,i,j,k,iblock,iparm)
+              v2c(2,l,i,j,k,iblock)=v2c_all(2,l,i,j,k,iblock,iparm)
+            enddo
+            do l=1,ntermd_2(i,j,k,iblock)
+              do m=1,ntermd_2(i,j,k,iblock)
+                v2s(l,m,i,j,k,iblock)=v2s_all(l,m,i,j,k,iblock,iparm)
+              enddo
+            enddo
+          enddo
+        enddo
+      enddo
+      enddo
+! Restore parameters of the cumulants
+      do i=-nloctyp,nloctyp
+        do j=1,2
+          b1(j,i)=b1_all(j,i,iparm)
+          b1tilde(j,i)=b1tilde_all(j,i,iparm)
+          b2(j,i)=b2_all(j,i,iparm)
+        enddo
+        do j=1,2
+          do k=1,2
+            cc(k,j,i)=cc_all(k,j,i,iparm)
+            ctilde(k,j,i)=ctilde_all(k,j,i,iparm)
+            dd(k,j,i)=dd_all(k,j,i,iparm)
+            dtilde(k,j,i)=dtilde_all(k,j,i,iparm)
+            ee(k,j,i)=ee_all(k,j,i,iparm)
+          enddo
+        enddo
+      enddo
+! Restore the parameters of electrostatic interactions
+      do i=1,2
+        do j=1,2
+          app(j,i)=app_all(j,i,iparm)
+          bpp(j,i)=bpp_all(j,i,iparm)
+          ael6(j,i)=ael6_all(j,i,iparm)
+          ael3(j,i)=ael3_all(j,i,iparm)
+        enddo
+      enddo
+! Restore sidechain parameters
+      do i=1,ntyp
+        do j=1,ntyp
+          aa(j,i)=aa_all(j,i,iparm)
+          bb(j,i)=bb_all(j,i,iparm)
+          r0(j,i)=r0_all(j,i,iparm)
+          sigma(j,i)=sigma_all(j,i,iparm)
+          chi(j,i)=chi_all(j,i,iparm)
+          augm(j,i)=augm_all(j,i,iparm)
+          eps(j,i)=eps_all(j,i,iparm)
+        enddo
+      enddo
+      do i=1,ntyp
+        chip(i)=chip_all(i,iparm)
+        alp(i)=alp_all(i,iparm)
+      enddo
+! Restore the SCp parameters
+      do i=1,ntyp
+        do j=1,2
+          aad(i,j)=aad_all(i,j,iparm)
+          bad(i,j)=bad_all(i,j,iparm)
+        enddo
+      enddo
+! Restore disulfide-bond parameters
+      ebr=ebr_all(iparm)
+      d0cm=d0cm_all(iparm)
+      akcm=akcm_all(iparm)
+      akth=akth_all(iparm)
+      akct=akct_all(iparm)
+      v1ss=v1ss_all(iparm)
+      v2ss=v2ss_all(iparm)
+      v3ss=v3ss_all(iparm)
+! Restore SC-backbone correlation parameters
+      do i=-nsccortyp,nsccortyp
+       do j=-nsccortyp,nsccortyp
+
+      nterm_sccor(j,i)=nterm_sccor_all(j,i,iparm)
+        do l=1,3
+           do k=1,nterm_sccor(j,i)
+            v1sccor(k,l,j,i)=v1sccor_all(k,l,j,i,iparm)
+            v2sccor(k,l,j,i)=v2sccor_all(k,l,j,i,iparm)
+           enddo
+          enddo
+        enddo
+      enddo
+      return
+      end subroutine restore_parm
+!--------------------------------------------------------------------------
+! make_ensemble1.F
+!--------------------------------------------------------------------------
+      subroutine make_ensembles(islice,*)
+! construct the conformational ensembles at REMD temperatures
+      use geometry_data, only:c
+      use io_base, only:ilen
+      use io_wham, only:pdboutW
+!      implicit none
+!      include "DIMENSIONS"
+!      include "DIMENSIONS.ZSCOPT"
+!      include "DIMENSIONS.FREE"
+#ifdef MPI
+      include "mpif.h"
+!      include "COMMON.MPI"
+      integer :: ierror,errcode,status(MPI_STATUS_SIZE) 
+#endif
+!      include "COMMON.IOUNITS"
+!      include "COMMON.CONTROL"
+!      include "COMMON.FREE"
+!      include "COMMON.ENERGIES"
+!      include "COMMON.FFIELD"
+!      include "COMMON.INTERACT"
+!      include "COMMON.SBRIDGE"
+!      include "COMMON.CHAIN"
+!      include "COMMON.PROTFILES"
+!      include "COMMON.PROT"
+      real(kind=4) :: csingle(3,nres*2)
+      real(kind=8),dimension(6) :: fT,fTprim,fTbis
+      real(kind=8) :: quot,quotl1,quotl,kfacl,&
+        eprim,ebis,temper,kfac=2.4d0,T0=300.0d0
+      real(kind=8) :: etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,&
+            escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,&
+            eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt
+      integer :: i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
+      real(kind=8) :: qfree,sumprob,eini,efree,rmsdev
+      character(len=80) :: bxname
+      character(len=2) :: licz1,licz2
+      character(len=3) :: licz3,licz4
+      character(len=5) :: ctemper
+!el      integer ilen
+!el      external ilen
+      real(kind=4) :: Fdimless(MaxStr),Fdimless_(MaxStr)
+      real(kind=8) :: enepot(MaxStr)
+      integer :: iperm(MaxStr)
+      integer :: islice
+      integer,dimension(0:nprocs) :: scount_
+#ifdef MPI
+      if (me.eq.Master) then
+#endif
+      write (licz2,'(bz,i2.2)') islice
+      if (nslice.eq.1) then
+        if (.not.separate_parset) then
+          bxname = prefix(:ilen(prefix))//".bx"
+        else
+          write (licz3,'(bz,i3.3)') myparm
+          bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx"
+        endif
+      else
+        if (.not.separate_parset) then
+          bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
+        else
+          write (licz3,'(bz,i3.3)') myparm
+          bxname = prefix(:ilen(prefix))//"par_"//licz3// &
+            "_slice_"//licz2//".bx"
+        endif
+      endif
+      open (ientout,file=bxname,status="unknown",&
+        form="unformatted",access="direct",recl=lenrec1)
+#ifdef MPI
+      endif
+#endif
+      do iparm=1,iparm
+        if (iparm.ne.iparmprint) exit
+        call restore_parm(iparm)
+        do ib=1,nT_h(iparm)
+#ifdef DEBUG
+          write (iout,*) "iparm",iparm," ib",ib
+#endif
+          temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
+!          quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
+!          quotl=1.0d0
+!          kfacl=1.0d0
+!          do l=1,5
+!            quotl1=quotl
+!            quotl=quotl*quot
+!            kfacl=kfacl*kfac
+!            fT(l)=kfacl/(kfacl-1.0d0+quotl)
+!          enddo
+!el old rescale weights
+!
+!            if (rescale_mode.eq.1) then
+!              quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
+#if defined(FUNCTH)
+              tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
+              ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
+#elif defined(FUNCT)
+              ft(6)=quot
+#else
+              ft(6)=1.0d0
+#endif
+!              quotl=1.0d0
+!              kfacl=1.0d0
+!              do l=1,5
+!                quotl1=quotl
+!                quotl=quotl*quot
+!                kfacl=kfacl*kfac
+!                fT(l)=kfacl/(kfacl-1.0d0+quotl)
+!              enddo
+!            else if (rescale_mode.eq.2) then
+!              quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
+!#if defined(FUNCTH)
+!              tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
+!              ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0
+!#elif defined(FUNCT)
+!              ft(6)=quot
+!#else 
+!              ft(6)=1.0d0
+!#endif
+!              quotl=1.0d0
+!              do l=1,5
+!                quotl=quotl*quot
+!                fT(l)=1.12692801104297249644d0/ &
+!                   dlog(dexp(quotl)+dexp(-quotl))
+!              enddo
+!              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
+!            else if (rescale_mode.eq.0) then
+!              do l=1,5
+!                fT(l)=0.0d0
+!              enddo
+!            else
+!              write (iout,*) &
+!              "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
+!              call flush(iout)
+!              return 1
+!            endif
+! el end old rescale weihgts
+          call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3))
+
+#ifdef MPI
+          do i=1,scount(me1)
+#else
+          do i=1,ntot(islice)
+#endif
+            evdw=enetb(1,i,iparm)
+!            evdw_t=enetb(21,i,iparm)
+            evdw_t=enetb(20,i,iparm)
+#ifdef SCP14
+!            evdw2_14=enetb(17,i,iparm)
+            evdw2_14=enetb(18,i,iparm)
+            evdw2=enetb(2,i,iparm)+evdw2_14
+#else
+            evdw2=enetb(2,i,iparm)
+            evdw2_14=0.0d0
+#endif
+#ifdef SPLITELE
+            ees=enetb(3,i,iparm)
+            evdw1=enetb(16,i,iparm)
+#else
+            ees=enetb(3,i,iparm)
+            evdw1=0.0d0
+#endif
+            ecorr=enetb(4,i,iparm)
+            ecorr5=enetb(5,i,iparm)
+            ecorr6=enetb(6,i,iparm)
+            eel_loc=enetb(7,i,iparm)
+            eello_turn3=enetb(8,i,iparm)
+            eello_turn4=enetb(9,i,iparm)
+            eello_turn6=enetb(10,i,iparm)
+            ebe=enetb(11,i,iparm)
+            escloc=enetb(12,i,iparm)
+            etors=enetb(13,i,iparm)
+            etors_d=enetb(14,i,iparm)
+            ehpb=enetb(15,i,iparm)
+
+            estr=enetb(17,i,iparm)
+!            estr=enetb(18,i,iparm)
+!            esccor=enetb(19,i,iparm)
+            esccor=enetb(21,i,iparm)
+!            edihcnstr=enetb(20,i,iparm)
+            edihcnstr=enetb(19,i,iparm)
+!#ifdef SPLITELE
+!            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
+!            +wvdwpp*evdw1 &
+!            +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
+!            +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
+!            +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
+!            +ft(2)*wturn3*eello_turn3 &
+!            +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc &
+!            +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+!            +wbond*estr
+!#else
+!            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
+!            +ft(1)*welec*(ees+evdw1) &
+!            +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
+!            +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
+!            +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
+!            +ft(2)*wturn3*eello_turn3 &
+!            +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr &
+!            +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+!            +wbond*estr
+!#endif
+
+#ifdef SPLITELE
+            etot=wsc*evdw+wscp*evdw2+welec*ees &
+            +wvdwpp*evdw1 &
+            +wang*ebe+wtor*etors+wscloc*escloc &
+            +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
+            +wcorr6*ecorr6+wturn4*eello_turn4 &
+            +wturn3*eello_turn3 &
+            +wturn6*eello_turn6+wel_loc*eel_loc &
+            +edihcnstr+wtor_d*etors_d+wsccor*esccor &
+            +wbond*estr
+#else
+            etot=wsc*evdw+wscp*evdw2 &
+            +welec*(ees+evdw1) &
+            +wang*ebe+wtor*etors+wscloc*escloc &
+            +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
+            +wcorr6*ecorr6+wturn4*eello_turn4 &
+            +wturn3*eello_turn3 &
+            +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
+            +wtor_d*etors_d+wsccor*esccor &
+            +wbond*estr
+#endif
+
+#ifdef MPI
+            Fdimless(i)= &
+              beta_h(ib,iparm)*etot-entfac(i)
+            potE(i,iparm)=etot
+#ifdef DEBUG
+            write (iout,*) i,indstart(me)+i-1,ib,&
+             1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),&
+             -entfac(i),Fdimless(i)
+#endif
+#else
+            Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
+            potE(i,iparm)=etot
+#endif
+          enddo   ! i
+#ifdef MPI
+          do i=1,scount(me1)
+            Fdimless_(i)=Fdimless(i)
+          enddo
+          call MPI_Gatherv(Fdimless_(1),scount(me),&
+           MPI_REAL,Fdimless(1),&
+           scount_(0),idispl(0),MPI_REAL,Master,&
+           WHAM_COMM, IERROR)
+#ifdef DEBUG
+          call MPI_Gatherv(potE(1,iparm),scount_(me),&
+           MPI_DOUBLE_PRECISION,potE(1,iparm),&
+           scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
+           WHAM_COMM, IERROR)
+          call MPI_Gatherv(entfac(1),scount(me),&
+           MPI_DOUBLE_PRECISION,entfac(1),&
+           scount_(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
+           WHAM_COMM, IERROR)
+#endif
+          if (me.eq.Master) then
+#ifdef DEBUG
+          write (iout,*) "The FDIMLESS array before sorting"
+          do i=1,ntot(islice)
+            write (iout,*) i,fdimless(i)
+          enddo
+#endif
+#endif
+          do i=1,ntot(islice)
+            iperm(i)=i
+          enddo
+          call mysort1(ntot(islice),Fdimless,iperm)
+#ifdef DEBUG
+          write (iout,*) "The FDIMLESS array after sorting"
+          do i=1,ntot(islice)
+            write (iout,*) i,iperm(i),fdimless(i)
+          enddo
+#endif
+          qfree=0.0d0
+          do i=1,ntot(islice)
+            qfree=qfree+exp(-fdimless(i)+fdimless(1))
+          enddo
+!          write (iout,*) "qfree",qfree
+          nlist=1
+          sumprob=0.0
+          do i=1,min0(ntot(islice),ensembles) 
+            sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree 
+#ifdef DEBUG
+            write (iout,*) i,ib,beta_h(ib,iparm),&
+             1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),&
+             potE(iperm(i),iparm),&
+             -entfac(iperm(i)),fdimless(i),sumprob
+#endif
+            if (sumprob.gt.0.99d0) goto 122
+            nlist=nlist+1
+          enddo  
+  122     continue
+#ifdef MPI
+          endif
+          call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,&
+             IERROR)
+          call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,&
+             IERROR)
+          do i=1,nlist
+            ii=iperm(i)
+            iproc=0
+            do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
+              iproc=iproc+1
+            enddo
+            if (iproc.ge.nprocs) then
+              write (iout,*) "Fatal error: processor out of range",iproc
+              call flush(iout)
+              if (bxfile) then
+                close (ientout)
+              else
+                close (ientout,status="delete")
+              endif
+              return 1
+            endif
+            ik=ii-indstart(iproc)+1
+            if (iproc.ne.Master) then
+              if (me.eq.iproc) then
+#ifdef DEBUG
+                write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,&
+                 " energy",potE(ik,iparm)
+#endif
+                call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,&
+                  Master,i,WHAM_COMM,IERROR)
+              else if (me.eq.Master) then
+                call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,&
+                  WHAM_COMM,STATUS,IERROR)
+              endif
+            else if (me.eq.Master) then
+              enepot(i)=potE(ik,iparm)
+            endif
+          enddo
+#else
+          do i=1,nlist
+            enepot(i)=potE(iperm(i),iparm)
+          enddo
+#endif
+#ifdef MPI
+          if (me.eq.Master) then
+#endif
+          write(licz3,'(bz,i3.3)') iparm
+          write(licz2,'(bz,i2.2)') islice
+          if (temper.lt.100.0d0) then
+            write(ctemper,'(f3.0)') temper
+          else if (temper.lt.1000.0) then
+            write (ctemper,'(f4.0)') temper
+          else
+            write (ctemper,'(f5.0)') temper
+          endif
+          if (nparmset.eq.1) then
+            if (separate_parset) then
+              write(licz4,'(bz,i3.3)') myparm
+              pdbname=prefix(:ilen(prefix))//"_par"//licz4
+            else
+              pdbname=prefix(:ilen(prefix))
+            endif
+          else
+            pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
+          endif
+          if (nslice.eq.1) then
+            pdbname=pdbname(:ilen(pdbname))//"_T_"// &
+              ctemper(:ilen(ctemper))//"pdb"
+          else
+            pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"// &
+              ctemper(:ilen(ctemper))//"pdb" 
+          endif
+          open(ipdb,file=pdbname)
+          do i=1,nlist
+            read (ientout,rec=iperm(i)) &
+              ((csingle(l,k),l=1,3),k=1,nres),&
+              ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
+              nss,(ihpb(k),jhpb(k),k=1,nss),&
+              eini,efree,rmsdev,iscor
+            do j=1,2*nres
+              do k=1,3
+                c(k,j)=csingle(k,j)
+              enddo
+            enddo
+            eini=fdimless(i)
+            call pdboutW(iperm(i),temper,eini,enepot(i),efree,rmsdev)
+          enddo
+#ifdef MPI
+          endif
+#endif
+        enddo     ! ib
+      enddo       ! iparm
+      if (bxfile) then
+        close(ientout)
+      else
+        close(ientout,status="delete")
+      endif
+      do i=0,nprocs
+        scount(i)=scount_(i)
+      enddo
+      return
+      end subroutine make_ensembles
+!--------------------------------------------------------------------------
+      subroutine mysort1(n, x, ipermut)
+!      implicit none
+      integer :: i,j,imax,ipm,n
+      real(kind=4) :: x(n)
+      integer :: ipermut(n)
+      real(kind=4) :: xtemp
+      do i=1,n
+        xtemp=x(i)
+        imax=i
+        do j=i+1,n
+          if (x(j).lt.xtemp) then
+            imax=j
+            xtemp=x(j)
+          endif
+        enddo
+        x(imax)=x(i)
+        x(i)=xtemp
+        ipm=ipermut(imax)
+        ipermut(imax)=ipermut(i)
+        ipermut(i)=ipm
+      enddo
+      return
+      end subroutine mysort1
+!--------------------------------------------------------------------------
+      subroutine alloc_enecalc_arrays(iparm)
+
+      use control_data
+      use geometry_data, only:maxlob
+      integer :: iparm
+!---------------------------
+! COMMON.ENERGIES form wham_data
+!      common /energies/
+      allocate(potE(MaxStr_Proc,iparm)) !(MaxStr_Proc,Max_Parm)
+      allocate(entfac(MaxStr_Proc)) !(MaxStr_Proc)
+      allocate(q(nQ+2,MaxStr_Proc)) !(MaxQ+2,MaxStr_Proc)
+      allocate(enetb(max_ene,MaxStr_Proc,iparm)) !(max_ene,MaxStr_Proc,Max_Parm)
+!
+! allocate ENECALC arrays
+!---------------------------
+! COMMON.ALLPARM
+!      common /allparm/
+      allocate(ww_all(max_eneW,iparm)) !(max_ene,max_parm) ! max_eneW
+      allocate(vbldp0_all(iparm),akp_all(nParmSet)) !(max_parm)
+      allocate(vbldsc0_all(maxbondterm,ntyp,iparm),&
+        aksc_all(maxbondterm,ntyp,iparm),&
+        abond0_all(maxbondterm,ntyp,iparm)) !(maxbondterm,ntyp,max_parm)
+      allocate(a0thet_all(-ntyp:ntyp,iparm)) !(-ntyp:ntyp,max_parm)
+      allocate(athet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm),&
+        bthet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm)) !(2,-ntyp:ntyp,-1:1,-1:1,max_parm)
+      allocate(polthet_all(0:3,-ntyp:ntyp,iparm)) !(0:3,-ntyp:ntyp,max_parm)
+      allocate(gthet_all(3,-ntyp:ntyp,iparm)) !(3,-ntyp:ntyp,max_parm)
+      allocate(theta0_all(-ntyp:ntyp,iparm),&
+        sig0_all(-ntyp:ntyp,iparm),sigc0_all(-ntyp:ntyp,nParmSet)) !(-ntyp:ntyp,max_parm)
+      allocate(aa0thet_all(-maxthetyp1:maxthetyp1,&
+        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,iparm))
+!(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
+      allocate(aathet_all(maxtheterm,-maxthetyp1:maxthetyp1,&
+        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,iparm))
+!(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
+      allocate(bbthet_all(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
+        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,iparm))
+      allocate(ccthet_all(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
+        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,iparm))
+      allocate(ddthet_all(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
+        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,iparm))
+      allocate(eethet_all(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
+        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,iparm))
+!(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,
+!     & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
+      allocate(ffthet_all1(maxdouble,maxdouble,maxtheterm3,&
+        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,&
+        -maxthetyp1:maxthetyp1,iparm))
+      allocate(ggthet_all1(maxdouble,maxdouble,maxtheterm3,&
+        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,&
+        -maxthetyp1:maxthetyp1,iparm))
+      allocate(ffthet_all2(maxdouble,maxdouble,maxtheterm3,&
+        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,&
+        -maxthetyp1:maxthetyp1,iparm))
+      allocate(ggthet_all2(maxdouble,maxdouble,maxtheterm3,&
+        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,&
+        -maxthetyp1:maxthetyp1,iparm))
+!(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
+!-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm)
+      allocate(dsc_all(ntyp1,iparm),dsc0_all(ntyp1,nParmSet)) !(ntyp1,max_parm)
+      allocate(bsc_all(maxlob,ntyp,iparm))
+!(maxlob,ntyp,max_parm)
+      allocate(censc_all(3,maxlob,-ntyp:ntyp,iparm)) !(3,maxlob,-ntyp:ntyp,max_parm)
+      allocate(gaussc_all(3,3,maxlob,-ntyp:ntyp,iparm)) !(3,3,maxlob,-ntyp:ntyp,max_parm)
+      allocate(sc_parmin_all(65,ntyp,iparm)) !(65,ntyp,max_parm)
+      allocate(v0_all(-maxtor:maxtor,-maxtor:maxtor,2,iparm))
+!(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
+      allocate(v1_all(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,iparm))
+      allocate(v2_all(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,iparm))
+!(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
+      allocate(vlor1_all(maxlor,maxtor,maxtor,iparm))
+      allocate(vlor2_all(maxlor,maxtor,maxtor,iparm))
+      allocate(vlor3_all(maxlor,maxtor,maxtor,iparm)) !(maxlor,maxtor,maxtor,max_parm)
+      allocate(v1c_all(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,&
+        -maxtor:maxtor,2,iparm))
+      allocate(v1s_all(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,&
+        -maxtor:maxtor,2,iparm))
+!(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
+      allocate(v2c_all(maxtermd_2,maxtermd_2,-maxtor:maxtor,&
+        -maxtor:maxtor,-maxtor:maxtor,2,iparm))
+!(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
+      allocate(v2s_all(maxtermd_2,maxtermd_2,-maxtor:maxtor,&
+        -maxtor:maxtor,-maxtor:maxtor,2,iparm))
+!(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
+      allocate(b1_all(2,-maxtor:maxtor,iparm))
+      allocate(b2_all(2,-maxtor:maxtor,iparm)) !(2,-maxtor:maxtor,max_parm)
+      allocate(cc_all(2,2,-maxtor:maxtor,iparm))
+      allocate(dd_all(2,2,-maxtor:maxtor,iparm))
+      allocate(ee_all(2,2,-maxtor:maxtor,iparm)) !(2,2,-maxtor:maxtor,max_parm)
+      allocate(ctilde_all(2,2,-maxtor:maxtor,iparm))
+      allocate(dtilde_all(2,2,-maxtor:maxtor,iparm)) !(2,2,-maxtor:maxtor,max_parm)
+      allocate(b1tilde_all(2,-maxtor:maxtor,iparm)) !(2,-maxtor:maxtor,max_parm)
+      allocate(app_all(2,2,iparm),bpp_all(2,2,nParmSet),&
+        ael6_all(2,2,iparm),ael3_all(2,2,nParmSet)) !(2,2,max_parm)
+      allocate(aad_all(ntyp,2,iparm),bad_all(ntyp,2,nParmSet)) !(ntyp,2,max_parm)
+      allocate(aa_all(ntyp,ntyp,iparm),bb_all(ntyp,ntyp,nParmSet),&
+        augm_all(ntyp,ntyp,iparm),eps_all(ntyp,ntyp,nParmSet),&
+        sigma_all(ntyp,ntyp,iparm),r0_all(ntyp,ntyp,nParmSet),&
+        chi_all(ntyp,ntyp,iparm)) !(ntyp,ntyp,max_parm)
+      allocate(chip_all(ntyp,iparm),alp_all(ntyp,nParmSet)) !(ntyp,max_parm)
+      allocate(ebr_all(iparm),d0cm_all(nParmSet),akcm_all(nParmSet),&
+        akth_all(iparm),akct_all(nParmSet),v1ss_all(nParmSet),&
+        v2ss_all(iparm),v3ss_all(nParmSet)) !(max_parm)
+      allocate(v1sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
+      allocate(v2sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
+!(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
+      allocate(nlob_all(ntyp1,iparm)) !(ntyp1,max_parm)
+      allocate(nlor_all(-maxtor:maxtor,-maxtor:maxtor,2,iparm))
+      allocate(nterm_all(-maxtor:maxtor,-maxtor:maxtor,2,iparm))
+!(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
+      allocate(ntermd1_all(-maxtor:maxtor,-maxtor:maxtor,&
+        -maxtor:maxtor,2,iparm))
+      allocate(ntermd2_all(-maxtor:maxtor,-maxtor:maxtor,&
+        -maxtor:maxtor,2,iparm))
+!(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
+      allocate(nbondterm_all(ntyp,iparm)) !(ntyp,max_parm)
+      allocate(ithetyp_all(-ntyp1:ntyp1,iparm)) !(-ntyp1:ntyp1,max_parm)
+      allocate(nthetyp_all(iparm),ntheterm_all(nParmSet),&
+        ntheterm2_all(iparm),ntheterm3_all(nParmSet),&
+        nsingle_all(iparm),&
+        ndouble_all(iparm),nntheterm_all(nParmSet)) !(max_parm)
+      allocate(nterm_sccor_all(-ntyp:ntyp,-ntyp:ntyp,iparm)) !(-ntyp:ntyp,-ntyp:ntyp,max_parm)
+!
+      end subroutine alloc_enecalc_arrays
+!--------------------------------------------------------------------------
+!--------------------------------------------------------------------------
+      end module ene_calc