source/wham/enecalc.F90

   1       module ene_calc
   2 !-----------------------------------------------------------------------------
   3       use io_units
   4       use wham_data
   5       use control_data, only: tor_mode
   6 !
   7       use geometry_data, only:nres,boxxsize,boxysize,boxzsize
   8       use energy_data
   9 !      use control_data, only:maxthetyp1
  10       use energy, only:etotal,enerprint,rescale_weights
  11 #ifdef MPI
  12       use MPI_data
  13 !      include "mpif.h"
  14 !      include "COMMON.MPI"
  15 #endif
  16       implicit none
  17 !-----------------------------------------------------------------------------
  18 ! COMMON.ALLPARM
  19 !      common /allparm/
  20       real(kind=8),dimension(:,:),allocatable :: ww_all !(max_ene,max_parm) ! max_eneW
  21       real(kind=8),dimension(:),allocatable :: vbldp0_all,akp_all !(max_parm)
  22       real(kind=8),dimension(:,:,:),allocatable :: vbldsc0_all,&
  23         aksc_all,abond0_all !(maxbondterm,ntyp,max_parm)
  24       real(kind=8),dimension(:,:),allocatable :: a0thet_all !(-ntyp:ntyp,max_parm)
  25       real(kind=8),dimension(:,:,:,:,:),allocatable :: athet_all,&
  26         bthet_all !(2,-ntyp:ntyp,-1:1,-1:1,max_parm)
  27       real(kind=8),dimension(:,:,:),allocatable :: polthet_all !(0:3,-ntyp:ntyp,max_parm)
  28       real(kind=8),dimension(:,:,:),allocatable :: gthet_all !(3,-ntyp:ntyp,max_parm)
  29       real(kind=8),dimension(:,:),allocatable :: theta0_all,&
  30         sig0_all,sigc0_all !(-ntyp:ntyp,max_parm)
  31       real(kind=8),dimension(:,:,:,:,:),allocatable :: aa0thet_all
  32 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
  33       real(kind=8),dimension(:,:,:,:,:,:),allocatable :: aathet_all
  34 !(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
  35       real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: bbthet_all,&
  36         ccthet_all,ddthet_all,eethet_all !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,
  37 !     & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
  38       real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: ffthet_all1,&
  39         ggthet_all1,ffthet_all2,ggthet_all2 !(maxdouble,maxdouble,maxtheterm3,
  40 !     &  -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm)
  41       real(kind=8),dimension(:,:),allocatable :: dsc_all,dsc0_all !(ntyp1,max_parm)
  42       real(kind=8),dimension(:,:,:),allocatable :: bsc_all !(maxlob,ntyp,max_parm)
  43       real(kind=8),dimension(:,:,:,:),allocatable :: censc_all !(3,maxlob,-ntyp:ntyp,max_parm)
  44       real(kind=8),dimension(:,:,:,:,:),allocatable :: gaussc_all !(3,3,maxlob,-ntyp:ntyp,max_parm)
  45       real(kind=8),dimension(:,:,:),allocatable :: sc_parmin_all !(65,ntyp,max_parm)
  46       real(kind=8),dimension(:,:,:,:),allocatable :: v0_all
  47 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
  48       real(kind=8),dimension(:,:,:,:,:),allocatable :: v1_all,&
  49         v2_all !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
  50       real(kind=8),dimension(:,:,:,:),allocatable :: vlor1_all,&
  51         vlor2_all,vlor3_all !(maxlor,maxtor,maxtor,max_parm)
  52       real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v1c_all,&
  53         v1s_all !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
  54       real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v2c_all
  55 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
  56       real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v2s_all
  57 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
  58       real(kind=8),dimension(:,:,:),allocatable :: b1_all,b2_all !(2,-maxtor:maxtor,max_parm)
  59       real(kind=8),dimension(:,:,:,:),allocatable :: cc_all,dd_all,&
  60         ee_all !(2,2,-maxtor:maxtor,max_parm)
  61       real(kind=8),dimension(:,:,:,:),allocatable :: ctilde_all,&
  62         dtilde_all !(2,2,-maxtor:maxtor,max_parm)
  63       real(kind=8),dimension(:,:,:),allocatable :: b1tilde_all !(2,-maxtor:maxtor,max_parm)
  64       real(kind=8),dimension(:,:,:),allocatable :: app_all,bpp_all,&
  65         ael6_all,ael3_all !(2,2,max_parm)
  66       real(kind=8),dimension(:,:,:),allocatable :: aad_all,&
  67         bad_all !(ntyp,2,max_parm)
  68       real(kind=8),dimension(:,:,:),allocatable :: aa_all,bb_all,&
  69         augm_all,eps_all,sigma_all,r0_all,chi_all !(ntyp,ntyp,max_parm)
  70       real(kind=8),dimension(:,:),allocatable :: chip_all,alp_all !(ntyp,max_parm)
  71       real(kind=8),dimension(:),allocatable :: ebr_all,d0cm_all,&
  72         akcm_all,akth_all,akct_all,v1ss_all,v2ss_all,v3ss_all !(max_parm)
  73       real(kind=8),dimension(:,:,:,:,:),allocatable :: v1sccor_all,&
  74         v2sccor_all !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
  75       integer,dimension(:,:),allocatable :: nlob_all !(ntyp1,max_parm)
  76       integer,dimension(:,:,:,:),allocatable :: nlor_all,&
  77         nterm_all !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
  78       integer,dimension(:,:,:,:,:),allocatable :: ntermd1_all,&
  79         ntermd2_all !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
  80       integer,dimension(:,:),allocatable :: nbondterm_all !(ntyp,max_parm)
  81       integer,dimension(:,:),allocatable :: ithetyp_all !(-ntyp1:ntyp1,max_parm)
  82       integer,dimension(:),allocatable :: nthetyp_all,ntheterm_all,&
  83         ntheterm2_all,ntheterm3_all,nsingle_all,ndouble_all,&
  84         nntheterm_all !(max_parm)
  85       integer,dimension(:,:,:),allocatable :: nterm_sccor_all !(-ntyp:ntyp,-ntyp:ntyp,max_parm)
  86 !-----------------------------------------------------------------------------
  87 !
  88 !
  89 !-----------------------------------------------------------------------------
  90       contains
  91 !-----------------------------------------------------------------------------
  92       subroutine enecalc(islice,*)
  93
  94       use names
  95       use control_data, only:indpdb
  96       use geometry_data, only:c,phi,theta,alph,omeg,deg2rad,anatemp,&
  97                               vbld,rad2deg,dc_norm,dc,vbld_inv
  98       use io_base, only:gyrate!,briefout
  99       use geometry, only:int_from_cart1,returnbox
 100       use io_wham, only:pdboutW
 101       use io_database, only:opentmp
 102       use conform_compar, only:qwolynes,rmsnat
 103 !      include "DIMENSIONS"
 104 !      include "DIMENSIONS.ZSCOPT"
 105 !      include "DIMENSIONS.FREE"
 106 #ifdef MPI
 107 !      use MPI_data
 108       include "mpif.h"
 109 !      include "COMMON.MPI"
 110 #endif
 111 !      include "COMMON.CHAIN"
 112 !      include "COMMON.IOUNITS"
 113 !      include "COMMON.PROTFILES"
 114 !      include "COMMON.NAMES"
 115 !      include "COMMON.VAR"
 116 !      include "COMMON.SBRIDGE"
 117 !      include "COMMON.GEO"
 118 !      include "COMMON.FFIELD"
 119 !      include "COMMON.ENEPS"
 120 !      include "COMMON.LOCAL"
 121 !      include "COMMON.WEIGHTS"
 122 !      include "COMMON.INTERACT"
 123 !      include "COMMON.FREE"
 124 !      include "COMMON.ENERGIES"
 125 !      include "COMMON.CONTROL"
 126 !      include "COMMON.TORCNSTR"
 127 !      implicit none
 128 #ifdef MPI
 129       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 130 #endif
 131       character(len=64) :: nazwa
 132       character(len=80) :: bxname
 133       character(len=3) :: liczba
 134 !el      real(kind=8) :: qwolynes
 135 !el      external qwolynes
 136       integer :: errmsg_count,maxerrmsg_count=100
 137 !el      real(kind=8) :: rmsnat,gyrate
 138 !el      external rmsnat,gyrate
 139       real(kind=8) :: tole=1.0d-1
 140       integer i,itj,ii,iii,j,k,l,licz
 141       integer ir,ib,ipar,iparm
 142       integer iscor,islice
 143       real(kind=4) :: csingle(3,nres*2)
 144       real(kind=8) :: energ
 145       real(kind=8) :: temp
 146 !el      integer ilen,iroof
 147 !el      external ilen,iroof
 148       real(kind=8) :: energia(0:n_ene),rmsdev,efree,eini
 149 !el      real(kind=8) :: energia(0:max_ene),rmsdev,efree,eini
 150       real(kind=8) :: fT(6),quot,quotl,kfacl,kfac=2.4d0,T0=3.0d2
 151       real(kind=8) :: tt
 152       integer :: snk_p(MaxR,MaxT_h,nParmSet)!Max_parm)
 153       logical :: lerr
 154       character(len=64) :: bprotfile_temp
 155
 156 !      integer :: rec
 157       integer,dimension(0:nprocs) :: scount_
 158 !el      real(kind=8) :: rmsnat
 159
 160       rescale_mode=rescale_modeW
 161
 162       call opentmp(islice,ientout,bprotfile_temp)
 163       iii=0
 164       ii=0
 165 !el
 166 !      iparm=1
 167       errmsg_count=0
 168       write (iout,*) "enecalc: nparmset ",nparmset
 169 #ifdef MPI
 170       do iparm=1,nParmSet
 171         do ib=1,nT_h(iparm)
 172           do i=1,nR(ib,iparm)
 173             snk_p(i,ib,iparm)=0
 174           enddo
 175         enddo
 176       enddo
 177       do i=indstart(me1),indend(me1)
 178 write(iout,*)"enecalc_ i indstart",i,indstart(me1),indend(me1)
 179 #else
 180       do iparm=1,nParmSet
 181         do ib=1,nT_h(iparm)
 182           do i=1,nR(ib,iparm)
 183             snk(i,ib,iparm)=0
 184           enddo
 185         enddo
 186       enddo
 187       do i=1,ntot
 188 write(iout,*)"enecalc_ i ntot",i,ntot
 189 #endif
 190         read(ientout,rec=i,err=101) &
 191           ((csingle(l,k),l=1,3),k=1,nres),&
 192           ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
 193           nss,(ihpb(k),jhpb(k),k=1,nss),&
 194           eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
 195 !el debug
 196 !write(iout,*)"co wczytuje"
 197 !          write(iout,*)((csingle(l,k),l=1,3),k=1,nres),&
 198 !          ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
 199 !          nss,(ihpb(k),jhpb(k),k=1,nss),&
 200 !          eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
 201 !el --------
 202
 203 !write(iout,*)"ipar",ib,ipar,1.0d0/(beta_h(ib,ipar)*1.987D-3)
 204          if (indpdb.gt.0) then
 205            do k=1,nres
 206              do l=1,3
 207                c(l,k)=csingle(l,k)
 208              enddo
 209            enddo
 210            do k=nnt,nct
 211              do l=1,3
 212                c(l,k+nres)=csingle(l,k+nres)
 213              enddo
 214            enddo
 215            call returnbox
 216            do k=1,nres
 217              do l=1,3
 218                csingle(l,k)=c(l,k)
 219              enddo
 220            enddo
 221            do k=nnt,nct
 222              do l=1,3
 223                csingle(l,k+nres)=c(l,k+nres)
 224              enddo
 225            enddo
 226            anatemp= 1.0d0/(beta_h(ib,ipar)*1.987D-3)
 227            q(nQ+1,iii+1)=rmsnat(iii+1)
 228          endif
 229          q(nQ+2,iii+1)=gyrate(iii+1)
 230 ! write(iout,*)"wczyt",anatemp,q(nQ+2,iii+1) !el
 231 !        fT=T0*beta_h(ib,ipar)*1.987D-3
 232 !        ft=2.0d0/(1.0d0+1.0d0/(T0*beta_h(ib,ipar)*1.987D-3))
 233 ! EL start old rescale
 234 !        if (rescale_modeW.eq.1) then
 235 !          quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
 236 !#if defined(FUNCTH)
 237 !          tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
 238 !          ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
 239 !#elif defined(FUNCT)
 240 !          ft(6)=quot
 241 !#else
 242 !          ft(6)=1.0d0
 243 !#endif
 244 !          quotl=1.0d0
 245 !          kfacl=1.0d0
 246 !          do l=1,5
 247 !            quotl=quotl*quot
 248 !            kfacl=kfacl*kfac
 249 !            fT(l)=kfacl/(kfacl-1.0d0+quotl)
 250 !          enddo
 251 !        else if (rescale_modeW.eq.2) then
 252 !          quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
 253 !#if defined(FUNCTH)
 254 !          tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
 255 !          ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
 256 !#elif defined(FUNCT)
 257 !          ft(6)=quot
 258 !#else
 259 !          ft(6)=1.0d0
 260 !#endif
 261 !          quotl=1.0d0
 262 !          do l=1,5
 263 !            quotl=quotl*quot
 264 !            fT(l)=1.12692801104297249644d0/ &
 265 !               dlog(dexp(quotl)+dexp(-quotl))
 266 !          enddo
 267 !        else if (rescale_modeW.eq.0) then
 268 !          do l=1,5
 269 !            fT(l)=1.0d0
 270 !          enddo
 271 !        else
 272 !          write (iout,*) "Error in ECECALC: wrong RESCALE_MODE",&
 273 !           rescale_modeW
 274 !          call flush(iout)
 275 !          return 1
 276 !        endif
 277 !EL end old rescele
 278 !        write (iout,*) "T",1.0d0/(beta_h(ib,ipar)*1.987D-3)," T0",T0,
 279 !     &   " kfac",kfac,"quot",quot," fT",fT
 280 #ifdef DEBUG
 281             write(iout,*)"weights"
 282             write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
 283             wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
 284             wtor_d,wsccor,wbond
 285 #endif
 286
 287         do j=1,2*nres
 288           do k=1,3
 289             c(k,j)=csingle(k,j)
 290           enddo
 291         enddo
 292         call int_from_cart1(.false.)
 293         ii=ii+1
 294
 295 !        call rescale_weights(1.0d0/(beta_h(ib,ipar)*1.987D-3))
 296        do iparm=1,nparmset
 297 #ifdef DEBUG
 298             write (iout,*) "before restore w=",1.0d0/(beta_h(ib,ipar)*1.987D-3)
 299             write(iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
 300             wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
 301             wtor_d,wsccor,wbond
 302 #endif
 303         call restore_parm(iparm)
 304         call rescale_weights(1.0d0/(beta_h(ib,ipar)*1.987D-3))
 305 #ifdef DEBUG
 306             write (iout,*) "before etot w=",1.0d0/(beta_h(ib,ipar)*1.987D-3)
 307             write(iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
 308             wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
 309             wtor_d,wsccor,wbond
 310 #endif
 311 !        call etotal(energia(0),fT)
 312         call etotal(energia(0))
 313 !write(iout,*)"check c and dc after etotal",1.0d0/(0.001987*beta_h(ib,ipar))
 314 !do k=1,2*nres+2
 315 !write(iout,*)k,"c=",(c(l,k),l=1,3)
 316 !write(iout,*)k,"dc=",(dc(l,k),l=1,3)
 317 !write(iout,*)k,"dc_norm=",(dc_norm(l,k),l=1,3)
 318 !enddo
 319 !do k=1,nres*2
 320 !write(iout,*)k,"vbld=",vbld(k)
 321 !write(iout,*)k,"vbld_inv=",vbld_inv(k)
 322 !enddo
 323
 324 !write(iout,*)"energia",(energia(j),j=0,n_ene)
 325 !write(iout,*)"enerprint tuz po call etotal"
 326         call enerprint(energia(0))
 327 #ifdef DEBUG
 328         write (iout,*) "Conformation",i
 329           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 330           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 331 !        call enerprint(energia(0),fT)
 332         call enerprint(energia(0))
 333         write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,49)
 334         write (iout,*) "ftors",ftors
 335 !el        call briefout(i,energia(0))
 336         temp=1.0d0/(beta_h(ib,ipar)*1.987D-3)
 337         write (iout,*) "temp", temp
 338         call pdboutW(i,temp,energia(0),energia(0),0.0d0,0.0d0)
 339 #endif
 340         if (energia(0).ge.1.0d20) then
 341           write (iout,*) "NaNs detected in some of the energy",&
 342            " components for conformation",ii+1
 343           write (iout,*) "The Cartesian geometry is:"
 344           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 345           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 346           write (iout,*) "The internal geometry is:"
 347 !          call intout
 348 !        call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
 349           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 350           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 351           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 352           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 353           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 354           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 355           write (iout,*) "The components of the energy are:"
 356 !          call enerprint(energia(0),fT)
 357           call enerprint(energia(0))
 358           write (iout,*) &
 359             "This conformation WILL NOT be added to the database."
 360           call flush(iout)
 361           goto 121
 362         else
 363 #ifdef DEBUG
 364           if (ipar.eq.iparm) write (iout,*) i,iparm,&
 365             1.0d0/(beta_h(ib,ipar)*1.987D-3),eini,energia(0)
 366 #endif
 367           if (ipar.eq.iparm .and. einicheck.gt.0 .and. &
 368             dabs(eini-energia(0)).gt.tole) then
 369             if (errmsg_count.le.maxerrmsg_count) then
 370               write (iout,'(2a,2e15.5,a,2i8,a,f8.1)') &
 371                "Warning: energy differs remarkably from ",&
 372                " the value read in: ",energia(0),eini," point",&
 373                iii+1,indstart(me1)+iii," T",&
 374                1.0d0/(1.987D-3*beta_h(ib,ipar))
 375 !              call intout
 376               call pdboutW(indstart(me1)+iii,&
 377        1.0d0/(1.987D-3*beta_h(ib,ipar)),&
 378        energia(0),eini,0.0d0,0.0d0)
 379           write (iout,*) "The Cartesian geometry is:"
 380           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 381           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 382           write (iout,*) "The internal geometry is:"
 383 !          call intout
 384 !        call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
 385           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 386           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 387           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 388           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 389           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 390           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 391               call enerprint(energia(0))
 392 !              call enerprint(energia(0),fT)
 393               errmsg_count=errmsg_count+1
 394               if (errmsg_count.gt.maxerrmsg_count) &
 395                 write (iout,*) "Too many warning messages"
 396               if (einicheck.gt.1) then
 397                 write (iout,*) "Calculation stopped."
 398                 call flush(iout)
 399 #ifdef MPI
 400                 call MPI_Abort(WHAM_COMM,IERROR,ERRCODE)
 401 #endif
 402                 call flush(iout)
 403                 return 1
 404               endif
 405             endif
 406           endif
 407           potE(iii+1,iparm)=energia(0)
 408           do k=1,49
 409             enetb(k,iii+1,iparm)=energia(k)
 410           enddo
 411 !           write (iout,'(2i5,21f8.2)') "debug",k,iii+1,(enetb(k,iii+1,iparm),k=1,21)
 412 !           write (iout,*) "debug",k,iii+1,(enetb(k,iii+1,iparm),k=1,21)
 413 #ifdef DEBUG
 414           write (iout,'(2i5,f10.1,3e15.5)') i,iii,&
 415            1.0d0/(beta_h(ib,ipar)*1.987D-3),energia(0),eini,efree
 416 !          call enerprint(energia(0),fT)
 417 #endif
 418 #ifdef DEBUG
 419           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 420           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 421           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 422           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 423           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 424           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 425           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 426           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 427           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 428           write (iout,'(8f10.5)') (q(k,iii+1),k=1,nQ)
 429           write (iout,'(f10.5,i10)') rmsdev,iscor
 430 !          call enerprint(energia(0),fT)
 431           call enerprint(energia(0))
 432         call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
 433 #endif
 434         endif
 435
 436         enddo ! iparm
 437
 438         iii=iii+1
 439         if (q(1,iii).le.0.0d0 .and. indpdb.gt.0) q(1,iii)=qwolynes(0,0)
 440         write (ientout,rec=iii) &
 441          ((csingle(l,k),l=1,3),k=1,nres),&
 442          ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
 443          nss,(ihpb(k),jhpb(k),k=1,nss),&
 444          potE(iii,ipar),efree,rmsdev,(q(k,iii),k=1,nQ),iR,ib,ipar
 445 !        write (iout,'(2i5,2e15.5)') ii,iii,potE(iii,ipar),efree
 446 #ifdef MPI
 447         if (separate_parset) then
 448           snk_p(iR,ib,1)=snk_p(iR,ib,1)+1
 449         else
 450           snk_p(iR,ib,ipar)=snk_p(iR,ib,ipar)+1
 451         endif
 452 !        write (iout,*) "iii",iii," iR",iR," ib",ib," ipar",ipar,
 453 !     &   " snk",snk_p(iR,ib,ipar)
 454 #else
 455         snk(iR,ib,ipar,islice)=snk(iR,ib,ipar,islice)+1
 456 #endif
 457   121   continue
 458       enddo
 459 #ifdef MPI
 460       scount(me)=iii
 461       write (iout,*) "Me",me," scount",scount(me)
 462       call flush(iout)
 463 !  Master gathers updated numbers of conformations written by all procs.
 464       call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount_(0), 1, &
 465         MPI_INTEGER, WHAM_COMM, IERROR)
 466       indstart(0)=1
 467       indend(0)=scount_(0)
 468       do i=1, Nprocs-1
 469         indstart(i)=indend(i-1)+1
 470         indend(i)=indstart(i)+scount_(i)-1
 471       enddo
 472       write (iout,*)
 473       write (iout,*) "Revised conformation counts"
 474       do i=0,nprocs1-1
 475         write (iout,'(a,i5,a,i7,a,i7,a,i7)') &
 476           "Processor",i," indstart",indstart(i),&
 477           " indend",indend(i)," count",scount_(i)
 478       enddo
 479       call flush(iout)
 480       call MPI_AllReduce(snk_p(1,1,1),snk(1,1,1,islice),&
 481         MaxR*MaxT_h*nParmSet,&
 482         MPI_INTEGER,MPI_SUM,WHAM_COMM,IERROR)
 483 #endif
 484       stot(islice)=0
 485       do iparm=1,nParmSet
 486         do ib=1,nT_h(iparm)
 487           do i=1,nR(ib,iparm)
 488             stot(islice)=stot(islice)+snk(i,ib,iparm,islice)
 489           enddo
 490         enddo
 491       enddo
 492       write (iout,*) "Revised SNK"
 493       do iparm=1,nParmSet
 494         do ib=1,nT_h(iparm)
 495           write (iout,'("Param",i3," Temp",f6.1,3x,32i8)') &
 496            iparm,1.0d0/(1.987D-3*beta_h(ib,iparm)),&
 497            (snk(i,ib,iparm,islice),i=1,nR(ib,iparm))
 498           write (iout,*) "snk_p",(snk_p(i,ib,iparm),i=1,nR(ib,iparm))
 499         enddo
 500       enddo
 501       write (iout,'("Total",i10)') stot(islice)
 502       call flush(iout)
 503       do i=0,nprocs
 504         scount(i)=scount_(i)
 505       enddo
 506       return
 507   101 write (iout,*) "Error in scratchfile."
 508       call flush(iout)
 509 !el#undef DEBUG
 510       return 1
 511       end subroutine enecalc
 512 !------------------------------------------------------------------------------
 513       logical function conf_check(ii,iprint)
 514
 515       use geometry_data, only:c,phi,theta,alph,omeg,deg2rad,rad2deg,vbld
 516       use geometry, only:int_from_cart1
 517 !      include "DIMENSIONS"
 518 !      include "DIMENSIONS.ZSCOPT"
 519 !      include "DIMENSIONS.FREE"
 520 #ifdef MPI
 521 !      use MPI_data
 522       include "mpif.h"
 523 !      include "COMMON.MPI"
 524 #endif
 525 !      include "COMMON.CHAIN"
 526 !      include "COMMON.IOUNITS"
 527 !      include "COMMON.PROTFILES"
 528 !      include "COMMON.NAMES"
 529 !      include "COMMON.VAR"
 530 !      include "COMMON.SBRIDGE"
 531 !      include "COMMON.GEO"
 532 !      include "COMMON.FFIELD"
 533 !      include "COMMON.ENEPS"
 534 !      include "COMMON.LOCAL"
 535 !      include "COMMON.WEIGHTS"
 536 !      include "COMMON.INTERACT"
 537 !      include "COMMON.FREE"
 538 !      include "COMMON.ENERGIES"
 539 !      include "COMMON.CONTROL"
 540 !      include "COMMON.TORCNSTR"
 541 !      implicit none
 542 #ifdef MPI
 543       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 544 #endif
 545       integer :: j,k,l,ii,itj,iprint,mnum,mnum1
 546       if (.not.check_conf) then
 547         conf_check=.true.
 548         return
 549       endif
 550       call int_from_cart1(.false.)
 551       do j=nnt+1,nct
 552         mnum=molnum(j)
 553         mnum1=molnum(j-1)
 554         if (mnum.ne.1) cycle
 555         if (itype(j-1,mnum1).ne.ntyp1_molec(mnum1) .and. itype(j,mnum).ne.ntyp1_molec(mnum) .and. &
 556           (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
 557           if (iprint.gt.0) &
 558           write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),&
 559             " for conformation",ii,mnum
 560           if (iprint.gt.1) then
 561             write (iout,*) "The Cartesian geometry is:"
 562             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 563             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 564             write (iout,*) "The internal geometry is:"
 565             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 566             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 567             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 568             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 569             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 570             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 571           endif
 572           if (iprint.gt.0) write (iout,*) &
 573             "This conformation WILL NOT be added to the database."
 574           conf_check=.false.
 575           return
 576         endif
 577       enddo
 578       do j=nnt,nct
 579         mnum=molnum(j)
 580        if (mnum.gt.3) cycle
 581         itj=itype(j,mnum)
 582         if (itype(j,mnum).ne.10 .and.itype(j,mnum).ne.ntyp1 .and. &
 583            (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then
 584           if (iprint.gt.0) &
 585           write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),&
 586            " for conformation",ii
 587           if (iprint.gt.1) then
 588             write (iout,*) "The Cartesian geometry is:"
 589             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 590             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 591             write (iout,*) "The internal geometry is:"
 592             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 593             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 594             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 595             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 596             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 597             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 598           endif
 599           if (iprint.gt.0) write (iout,*) &
 600             "This conformation WILL NOT be added to the database."
 601           conf_check=.false.
 602           return
 603         endif
 604       enddo
 605       do j=3,nres
 606         if (theta(j).le.0.0d0) then
 607           if (iprint.gt.0) &
 608           write (iout,*) "Zero theta angle(s) in conformation",ii
 609           if (iprint.gt.1) then
 610             write (iout,*) "The Cartesian geometry is:"
 611             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 612             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 613             write (iout,*) "The internal geometry is:"
 614             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 615             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 616             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 617             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 618             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 619             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 620           endif
 621           if (iprint.gt.0) write (iout,*) &
 622             "This conformation WILL NOT be added to the database."
 623           conf_check=.false.
 624           return
 625         endif
 626         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
 627       enddo
 628       conf_check=.true.
 629 !      write (iout,*) "conf_check passed",ii
 630       return
 631       end function conf_check
 632 !-----------------------------------------------------------------------------
 633 ! store_parm.F
 634 !-----------------------------------------------------------------------------
 635       subroutine store_parm(iparm)
 636 !
 637 ! Store parameters of set IPARM
 638 ! valence angles and the side chains and energy parameters.
 639 !
 640 !      implicit none
 641 !      include 'DIMENSIONS'
 642 !      include 'DIMENSIONS.ZSCOPT'
 643 !      include 'DIMENSIONS.FREE'
 644 !      include 'COMMON.IOUNITS'
 645 !      include 'COMMON.CHAIN'
 646 !      include 'COMMON.INTERACT'
 647 !      include 'COMMON.GEO'
 648 !      include 'COMMON.LOCAL'
 649 !      include 'COMMON.TORSION'
 650 !      include 'COMMON.FFIELD'
 651 !      include 'COMMON.NAMES'
 652 !      include 'COMMON.SBRIDGE'
 653 !      include 'COMMON.SCROT'
 654 !      include 'COMMON.SCCOR'
 655 !      include 'COMMON.ALLPARM'
 656       integer :: i,j,k,l,m,mm,iparm,ichir1,ichir2,iblock,iii
 657
 658       call alloc_enecalc_arrays(iparm)
 659 !el      allocate(ww_all(n_ene,iparm))
 660 ! Store weights
 661       ww_all(1,iparm)=wsc
 662       ww_all(2,iparm)=wscp
 663       ww_all(3,iparm)=welec
 664       ww_all(4,iparm)=wcorr
 665       ww_all(5,iparm)=wcorr5
 666       ww_all(6,iparm)=wcorr6
 667       ww_all(7,iparm)=wel_loc
 668       ww_all(8,iparm)=wturn3
 669       ww_all(9,iparm)=wturn4
 670       ww_all(10,iparm)=wturn6
 671       ww_all(11,iparm)=wang
 672       ww_all(12,iparm)=wscloc
 673       ww_all(13,iparm)=wtor
 674       ww_all(14,iparm)=wtor_d
 675       ww_all(15,iparm)=wstrain
 676       ww_all(16,iparm)=wvdwpp
 677       ww_all(17,iparm)=wbond
 678       ww_all(19,iparm)=wsccor
 679       ww_all(42,iparm)=wcatprot
 680       ww_all(41,iparm)=wcatcat
 681 ! Store bond parameters
 682       vbldp0_all(iparm)=vbldp0
 683       akp_all(iparm)=akp
 684       do i=1,ntyp
 685         nbondterm_all(i,iparm)=nbondterm(i)
 686         do j=1,nbondterm(i)
 687           vbldsc0_all(j,i,iparm)=vbldsc0(j,i)
 688           aksc_all(j,i,iparm)=aksc(j,i)
 689           abond0_all(j,i,iparm)=abond0(j,i)
 690         enddo
 691       enddo
 692 ! Store bond angle parameters
 693 #ifdef CRYST_THETA
 694       do i=-ntyp,ntyp
 695         a0thet_all(i,iparm)=a0thet(i)
 696         do ichir1=-1,1
 697         do ichir2=-1,1
 698         do j=1,2
 699           athet_all(j,i,ichir1,ichir2,iparm)=athet(j,i,ichir1,ichir2)
 700           bthet_all(j,i,ichir1,ichir2,iparm)=bthet(j,i,ichir1,ichir2)
 701         enddo
 702         enddo
 703         enddo
 704         do j=0,3
 705           polthet_all(j,i,iparm)=polthet(j,i)
 706         enddo
 707         do j=1,3
 708           gthet_all(j,i,iparm)=gthet(j,i)
 709         enddo
 710         theta0_all(i,iparm)=theta0(i)
 711         sig0_all(i,iparm)=sig0(i)
 712         sigc0_all(i,iparm)=sigc0(i)
 713       enddo
 714 #else
 715       IF (tor_mode.eq.0) THEN
 716       nthetyp_all(iparm)=nthetyp
 717       ntheterm_all(iparm)=ntheterm
 718       ntheterm2_all(iparm)=ntheterm2
 719       ntheterm3_all(iparm)=ntheterm3
 720       nsingle_all(iparm)=nsingle
 721       ndouble_all(iparm)=ndouble
 722       nntheterm_all(iparm)=nntheterm
 723       do i=-ntyp,ntyp
 724         ithetyp_all(i,iparm)=ithetyp(i)
 725       enddo
 726       do iblock=1,2
 727       do i=-nthetyp-1,nthetyp+1
 728         do j=-nthetyp-1,nthetyp+1
 729           do k=-nthetyp-1,nthetyp+1
 730             aa0thet_all(i,j,k,iblock,iparm)=aa0thet(i,j,k,iblock)
 731             do l=1,ntheterm
 732               aathet_all(l,i,j,k,iblock,iparm)=aathet(l,i,j,k,iblock)
 733             enddo
 734             do l=1,ntheterm2
 735               do m=1,nsingle
 736                 bbthet_all(m,l,i,j,k,iblock,iparm)= &
 737       bbthet(m,l,i,j,k,iblock)
 738                 ccthet_all(m,l,i,j,k,iblock,iparm)= &
 739       ccthet(m,l,i,j,k,iblock)
 740                 ddthet_all(m,l,i,j,k,iblock,iparm)= &
 741       ddthet(m,l,i,j,k,iblock)
 742                 eethet_all(m,l,i,j,k,iblock,iparm)= &
 743       eethet(m,l,i,j,k,iblock)
 744               enddo
 745             enddo
 746             do l=1,ntheterm3
 747               do m=1,ndouble
 748                 do mm=1,ndouble
 749                 if (iblock.eq.1) then
 750                  ffthet_all1(mm,m,l,i,j,k,iparm)=&
 751          ffthet(mm,m,l,i,j,k,iblock)
 752                  ggthet_all1(mm,m,l,i,j,k,iparm)=&
 753       ggthet(mm,m,l,i,j,k,iblock)
 754                   else
 755                  ffthet_all2(mm,m,l,i,j,k,iparm)=&
 756          ffthet(mm,m,l,i,j,k,iblock)
 757                  ggthet_all2(mm,m,l,i,j,k,iparm)=&
 758       ggthet(mm,m,l,i,j,k,iblock)
 759                   endif
 760                 enddo
 761               enddo
 762             enddo
 763           enddo
 764         enddo
 765       enddo
 766       enddo
 767       ELSE
 768       write(iout,*) "Need storing parameters"
 769       ENDIF
 770 #endif
 771 #ifdef CRYST_SC
 772 ! Store the sidechain rotamer parameters
 773       do i=-ntyp,ntyp
 774        iii=iabs(i)
 775 !!       write (iout,*) i,"storeparm1"
 776        if (i.eq.0) cycle
 777         nlob_all(iii,iparm)=nlob(iii)
 778         do j=1,nlob(iii)
 779           bsc_all(j,iii,iparm)=bsc(j,iii)
 780           do k=1,3
 781             censc_all(k,j,i,iparm)=censc(k,j,i)
 782           enddo
 783           do k=1,3
 784             do l=1,3
 785               gaussc_all(l,k,j,i,iparm)=gaussc(l,k,j,i)
 786             enddo
 787           enddo
 788         enddo
 789       enddo
 790 #else
 791       do i=1,ntyp
 792         do j=1,65
 793           sc_parmin_all(j,i,iparm)=sc_parmin(j,i)
 794         enddo
 795       enddo
 796 #endif
 797       IF (tor_mode.eq.0) THEN
 798 ! Store the torsional parameters
 799       do iblock=1,2
 800       do i=-ntortyp+1,ntortyp-1
 801         do j=-ntortyp+1,ntortyp-1
 802           v0_all(i,j,iblock,iparm)=v0(i,j,iblock)
 803           nterm_all(i,j,iblock,iparm)=nterm(i,j,iblock)
 804           nlor_all(i,j,iblock,iparm)=nlor(i,j,iblock)
 805           do k=1,nterm(i,j,iblock)
 806             v1_all(k,i,j,iblock,iparm)=v1(k,i,j,iblock)
 807             v2_all(k,i,j,iblock,iparm)=v2(k,i,j,iblock)
 808           enddo
 809           do k=1,nlor(i,j,iblock)
 810             vlor1_all(k,i,j,iparm)=vlor1(k,i,j)
 811             vlor2_all(k,i,j,iparm)=vlor2(k,i,j)
 812             vlor3_all(k,i,j,iparm)=vlor3(k,i,j)
 813           enddo
 814         enddo
 815       enddo
 816       enddo
 817 ! Store the double torsional parameters
 818       do iblock=1,2
 819       do i=-ntortyp+1,ntortyp-1
 820         do j=-ntortyp+1,ntortyp-1
 821           do k=-ntortyp+1,ntortyp-1
 822             ntermd1_all(i,j,k,iblock,iparm)=ntermd_1(i,j,k,iblock)
 823             ntermd2_all(i,j,k,iblock,iparm)=ntermd_2(i,j,k,iblock)
 824             do l=1,ntermd_1(i,j,k,iblock)
 825               v1c_all(1,l,i,j,k,iblock,iparm)=v1c(1,l,i,j,k,iblock)
 826               v1c_all(2,l,i,j,k,iblock,iparm)=v1c(2,l,i,j,k,iblock)
 827               v2c_all(1,l,i,j,k,iblock,iparm)=v2c(1,l,i,j,k,iblock)
 828               v2c_all(2,l,i,j,k,iblock,iparm)=v2c(2,l,i,j,k,iblock)
 829             enddo
 830             do l=1,ntermd_2(i,j,k,iblock)
 831               do m=1,ntermd_2(i,j,k,iblock)
 832                 v2s_all(l,m,i,j,k,iblock,iparm)=v2s(l,m,i,j,k,iblock)
 833               enddo
 834             enddo
 835           enddo
 836         enddo
 837       enddo
 838       enddo
 839 ! Store parameters of the cumulants
 840       do i=1,nloctyp
 841         do j=1,2
 842           b1_all(j,i,iparm)=b1(j,i)
 843           b1tilde_all(j,i,iparm)=b1tilde(j,i)
 844           b2_all(j,i,iparm)=b2(j,i)
 845         enddo
 846         do j=1,2
 847           do k=1,2
 848             cc_all(k,j,i,iparm)=cc(k,j,i)
 849             ctilde_all(k,j,i,iparm)=ctilde(k,j,i)
 850             dd_all(k,j,i,iparm)=dd(k,j,i)
 851             dtilde_all(k,j,i,iparm)=dtilde(k,j,i)
 852             ee_all(k,j,i,iparm)=ee(k,j,i)
 853           enddo
 854         enddo
 855       enddo
 856       ELSE
 857        write(iout,*) "NEED storing parameters"
 858       ENDIF
 859 ! Store the parameters of electrostatic interactions
 860       do i=1,2
 861         do j=1,2
 862           app_all(j,i,iparm)=app(j,i)
 863           bpp_all(j,i,iparm)=bpp(j,i)
 864           ael6_all(j,i,iparm)=ael6(j,i)
 865           ael3_all(j,i,iparm)=ael3(j,i)
 866         enddo
 867       enddo
 868 ! Store sidechain parameters
 869       do i=1,ntyp
 870         do j=1,ntyp
 871           aa_all(j,i,iparm)=aa_aq(j,i)
 872           bb_all(j,i,iparm)=bb_aq(j,i)
 873           r0_all(j,i,iparm)=r0(j,i)
 874           sigma_all(j,i,iparm)=sigma(j,i)
 875           chi_all(j,i,iparm)=chi(j,i)
 876           augm_all(j,i,iparm)=augm(j,i)
 877           eps_all(j,i,iparm)=eps(j,i)
 878         enddo
 879       enddo
 880       do i=1,ntyp
 881         chip_all(i,iparm)=chip(i)
 882         alp_all(i,iparm)=alp(i)
 883       enddo
 884 ! Store the SCp parameters
 885       do i=1,ntyp
 886         do j=1,2
 887           aad_all(i,j,iparm)=aad(i,j)
 888           bad_all(i,j,iparm)=bad(i,j)
 889         enddo
 890       enddo
 891 ! Store disulfide-bond parameters
 892       ebr_all(iparm)=ebr
 893       d0cm_all(iparm)=d0cm
 894       akcm_all(iparm)=akcm
 895       akth_all(iparm)=akth
 896       akct_all(iparm)=akct
 897       v1ss_all(iparm)=v1ss
 898       v2ss_all(iparm)=v2ss
 899       v3ss_all(iparm)=v3ss
 900 ! Store SC-backbone correlation parameters
 901       do i=-nsccortyp,nsccortyp
 902        do j=-nsccortyp,nsccortyp
 903
 904       nterm_sccor_all(j,i,iparm)=nterm_sccor(j,i)
 905 !      do i=1,20
 906 !        do j=1,20
 907          do l=1,3
 908           do k=1,nterm_sccor(j,i)
 909             v1sccor_all(k,l,j,i,iparm)=v1sccor(k,l,j,i)
 910             v2sccor_all(k,l,j,i,iparm)=v2sccor(k,l,j,i)
 911            enddo
 912           enddo
 913         enddo
 914       enddo
 915 write(iout,*)"end of store_parm"
 916       return
 917       end subroutine store_parm
 918 !--------------------------------------------------------------------------
 919       subroutine restore_parm(iparm)
 920 !
 921 ! Store parameters of set IPARM
 922 ! valence angles and the side chains and energy parameters.
 923 !
 924 !      implicit none
 925 !      include 'DIMENSIONS'
 926 !      include 'DIMENSIONS.ZSCOPT'
 927 !      include 'DIMENSIONS.FREE'
 928 !      include 'COMMON.IOUNITS'
 929 !      include 'COMMON.CHAIN'
 930 !      include 'COMMON.INTERACT'
 931 !      include 'COMMON.GEO'
 932 !      include 'COMMON.LOCAL'
 933 !      include 'COMMON.TORSION'
 934 !      include 'COMMON.FFIELD'
 935 !      include 'COMMON.NAMES'
 936 !      include 'COMMON.SBRIDGE'
 937 !      include 'COMMON.SCROT'
 938 !      include 'COMMON.SCCOR'
 939 !      include 'COMMON.ALLPARM'
 940       integer :: i,j,k,l,m,mm,iparm,ichir1,ichir2,iblock,iii
 941
 942 ! Restore weights
 943       wsc=ww_all(1,iparm)
 944       wscp=ww_all(2,iparm)
 945       welec=ww_all(3,iparm)
 946       wcorr=ww_all(4,iparm)
 947       wcorr5=ww_all(5,iparm)
 948       wcorr6=ww_all(6,iparm)
 949       wel_loc=ww_all(7,iparm)
 950       wturn3=ww_all(8,iparm)
 951       wturn4=ww_all(9,iparm)
 952       wturn6=ww_all(10,iparm)
 953       wang=ww_all(11,iparm)
 954       wscloc=ww_all(12,iparm)
 955       wtor=ww_all(13,iparm)
 956       wtor_d=ww_all(14,iparm)
 957       wstrain=ww_all(15,iparm)
 958       wvdwpp=ww_all(16,iparm)
 959       wbond=ww_all(17,iparm)
 960       wsccor=ww_all(19,iparm)
 961       wcatprot=ww_all(42,iparm)
 962       wcatcat=ww_all(41,iparm)
 963 ! Restore bond parameters
 964       vbldp0=vbldp0_all(iparm)
 965       akp=akp_all(iparm)
 966       do i=1,ntyp
 967         nbondterm(i)=nbondterm_all(i,iparm)
 968         do j=1,nbondterm(i)
 969           vbldsc0(j,i)=vbldsc0_all(j,i,iparm)
 970           aksc(j,i)=aksc_all(j,i,iparm)
 971           abond0(j,i)=abond0_all(j,i,iparm)
 972         enddo
 973       enddo
 974 ! Restore bond angle parameters
 975 #ifdef CRYST_THETA
 976       do i=-ntyp,ntyp
 977         a0thet(i)=a0thet_all(i,iparm)
 978         do ichir1=-1,1
 979         do ichir2=-1,1
 980         do j=1,2
 981           athet(j,i,ichir1,ichir2)=athet_all(j,i,ichir1,ichir2,iparm)
 982           bthet(j,i,ichir1,ichir2)=bthet_all(j,i,ichir1,ichir2,iparm)
 983         enddo
 984         enddo
 985         enddo
 986         do j=0,3
 987           polthet(j,i)=polthet_all(j,i,iparm)
 988         enddo
 989         do j=1,3
 990           gthet(j,i)=gthet_all(j,i,iparm)
 991         enddo
 992         theta0(i)=theta0_all(i,iparm)
 993         sig0(i)=sig0_all(i,iparm)
 994         sigc0(i)=sigc0_all(i,iparm)
 995       enddo
 996 #else
 997       IF (tor_mode.eq.0) THEN
 998       nthetyp=nthetyp_all(iparm)
 999       ntheterm=ntheterm_all(iparm)
1000       ntheterm2=ntheterm2_all(iparm)
1001       ntheterm3=ntheterm3_all(iparm)
1002       nsingle=nsingle_all(iparm)
1003       ndouble=ndouble_all(iparm)
1004       nntheterm=nntheterm_all(iparm)
1005       do i=-ntyp,ntyp
1006         ithetyp(i)=ithetyp_all(i,iparm)
1007       enddo
1008       do iblock=1,2
1009       do i=-nthetyp-1,nthetyp+1
1010         do j=-nthetyp-1,nthetyp+1
1011           do k=-nthetyp-1,nthetyp+1
1012             aa0thet(i,j,k,iblock)=aa0thet_all(i,j,k,iblock,iparm)
1013             do l=1,ntheterm
1014               aathet(l,i,j,k,iblock)=aathet_all(l,i,j,k,iblock,iparm)
1015             enddo
1016             do l=1,ntheterm2
1017               do m=1,nsingle
1018                 bbthet(m,l,i,j,k,iblock)= &
1019       bbthet_all(m,l,i,j,k,iblock,iparm)
1020                 ccthet(m,l,i,j,k,iblock)= &
1021       ccthet_all(m,l,i,j,k,iblock,iparm)
1022                 ddthet(m,l,i,j,k,iblock)= &
1023       ddthet_all(m,l,i,j,k,iblock,iparm)
1024                 eethet(m,l,i,j,k,iblock)= &
1025       eethet_all(m,l,i,j,k,iblock,iparm)
1026               enddo
1027             enddo
1028             do l=1,ntheterm3
1029               do m=1,ndouble
1030                 do mm=1,ndouble
1031                 if (iblock.eq.1) then
1032                  ffthet(mm,m,l,i,j,k,iblock)= &
1033       ffthet_all1(mm,m,l,i,j,k,iparm)
1034                  ggthet(mm,m,l,i,j,k,iblock)= &
1035       ggthet_all1(mm,m,l,i,j,k,iparm)
1036                 else
1037                  ffthet(mm,m,l,i,j,k,iblock)= &
1038       ffthet_all2(mm,m,l,i,j,k,iparm)
1039                  ggthet(mm,m,l,i,j,k,iblock)= &
1040       ggthet_all2(mm,m,l,i,j,k,iparm)
1041                 endif
1042                 enddo
1043               enddo
1044             enddo
1045           enddo
1046         enddo
1047       enddo
1048       enddo
1049       ELSE
1050       write (iout,*) "Need storing parameters"
1051       ENDIF
1052 #endif
1053 ! Restore the sidechain rotamer parameters
1054 #ifdef CRYST_SC
1055       do i=-ntyp,ntyp
1056         if (i.eq.0) cycle
1057         iii=iabs(i)
1058         nlob(iii)=nlob_all(iii,iparm)
1059         do j=1,nlob(iii)
1060           bsc(j,iii)=bsc_all(j,iii,iparm)
1061           do k=1,3
1062             censc(k,j,i)=censc_all(k,j,i,iparm)
1063           enddo
1064           do k=1,3
1065             do l=1,3
1066               gaussc(l,k,j,i)=gaussc_all(l,k,j,i,iparm)
1067             enddo
1068           enddo
1069         enddo
1070       enddo
1071 #else
1072       do i=1,ntyp
1073         do j=1,65
1074           sc_parmin(j,i)=sc_parmin_all(j,i,iparm)
1075         enddo
1076       enddo
1077 #endif
1078       IF (tor_mode.eq.0) THEN
1079 ! Restore the torsional parameters
1080       do iblock=1,2
1081       do i=-ntortyp+1,ntortyp-1
1082         do j=-ntortyp+1,ntortyp-1
1083           v0(i,j,iblock)=v0_all(i,j,iblock,iparm)
1084           nterm(i,j,iblock)=nterm_all(i,j,iblock,iparm)
1085           nlor(i,j,iblock)=nlor_all(i,j,iblock,iparm)
1086           do k=1,nterm(i,j,iblock)
1087             v1(k,i,j,iblock)=v1_all(k,i,j,iblock,iparm)
1088             v2(k,i,j,iblock)=v2_all(k,i,j,iblock,iparm)
1089           enddo
1090           do k=1,nlor(i,j,iblock)
1091             vlor1(k,i,j)=vlor1_all(k,i,j,iparm)
1092             vlor2(k,i,j)=vlor2_all(k,i,j,iparm)
1093             vlor3(k,i,j)=vlor3_all(k,i,j,iparm)
1094           enddo
1095         enddo
1096       enddo
1097       enddo
1098 ! Restore the double torsional parameters
1099       do iblock=1,2
1100       do i=-ntortyp+1,ntortyp-1
1101         do j=-ntortyp+1,ntortyp-1
1102           do k=-ntortyp+1,ntortyp-1
1103             ntermd_1(i,j,k,iblock)=ntermd1_all(i,j,k,iblock,iparm)
1104             ntermd_2(i,j,k,iblock)=ntermd2_all(i,j,k,iblock,iparm)
1105             do l=1,ntermd_1(i,j,k,iblock)
1106               v1c(1,l,i,j,k,iblock)=v1c_all(1,l,i,j,k,iblock,iparm)
1107               v1c(2,l,i,j,k,iblock)=v1c_all(2,l,i,j,k,iblock,iparm)
1108               v2c(1,l,i,j,k,iblock)=v2c_all(1,l,i,j,k,iblock,iparm)
1109               v2c(2,l,i,j,k,iblock)=v2c_all(2,l,i,j,k,iblock,iparm)
1110             enddo
1111             do l=1,ntermd_2(i,j,k,iblock)
1112               do m=1,ntermd_2(i,j,k,iblock)
1113                 v2s(l,m,i,j,k,iblock)=v2s_all(l,m,i,j,k,iblock,iparm)
1114               enddo
1115             enddo
1116           enddo
1117         enddo
1118       enddo
1119       enddo
1120
1121 ! Restore parameters of the cumulants
1122       do i=1,nloctyp
1123         do j=1,2
1124           b1(j,i)=b1_all(j,i,iparm)
1125           b1tilde(j,i)=b1tilde_all(j,i,iparm)
1126           b2(j,i)=b2_all(j,i,iparm)
1127         enddo
1128         do j=1,2
1129           do k=1,2
1130             cc(k,j,i)=cc_all(k,j,i,iparm)
1131             ctilde(k,j,i)=ctilde_all(k,j,i,iparm)
1132             dd(k,j,i)=dd_all(k,j,i,iparm)
1133             dtilde(k,j,i)=dtilde_all(k,j,i,iparm)
1134             ee(k,j,i)=ee_all(k,j,i,iparm)
1135           enddo
1136         enddo
1137       enddo
1138       ELSE
1139       write(iout,*) "need storing parameters"
1140       ENDIF
1141 ! Restore the parameters of electrostatic interactions
1142       do i=1,2
1143         do j=1,2
1144           app(j,i)=app_all(j,i,iparm)
1145           bpp(j,i)=bpp_all(j,i,iparm)
1146           ael6(j,i)=ael6_all(j,i,iparm)
1147           ael3(j,i)=ael3_all(j,i,iparm)
1148         enddo
1149       enddo
1150 ! Restore sidechain parameters
1151       do i=1,ntyp
1152         do j=1,ntyp
1153           aa_aq(j,i)=aa_all(j,i,iparm)
1154           bb_aq(j,i)=bb_all(j,i,iparm)
1155           r0(j,i)=r0_all(j,i,iparm)
1156           sigma(j,i)=sigma_all(j,i,iparm)
1157           chi(j,i)=chi_all(j,i,iparm)
1158           augm(j,i)=augm_all(j,i,iparm)
1159           eps(j,i)=eps_all(j,i,iparm)
1160         enddo
1161       enddo
1162       do i=1,ntyp
1163         chip(i)=chip_all(i,iparm)
1164         alp(i)=alp_all(i,iparm)
1165       enddo
1166 ! Restore the SCp parameters
1167       do i=1,ntyp
1168         do j=1,2
1169           aad(i,j)=aad_all(i,j,iparm)
1170           bad(i,j)=bad_all(i,j,iparm)
1171         enddo
1172       enddo
1173 ! Restore disulfide-bond parameters
1174       ebr=ebr_all(iparm)
1175       d0cm=d0cm_all(iparm)
1176       akcm=akcm_all(iparm)
1177       akth=akth_all(iparm)
1178       akct=akct_all(iparm)
1179       v1ss=v1ss_all(iparm)
1180       v2ss=v2ss_all(iparm)
1181       v3ss=v3ss_all(iparm)
1182 ! Restore SC-backbone correlation parameters
1183       do i=-nsccortyp,nsccortyp
1184        do j=-nsccortyp,nsccortyp
1185
1186       nterm_sccor(j,i)=nterm_sccor_all(j,i,iparm)
1187         do l=1,3
1188            do k=1,nterm_sccor(j,i)
1189             v1sccor(k,l,j,i)=v1sccor_all(k,l,j,i,iparm)
1190             v2sccor(k,l,j,i)=v2sccor_all(k,l,j,i,iparm)
1191            enddo
1192           enddo
1193         enddo
1194       enddo
1195       return
1196       end subroutine restore_parm
1197 !--------------------------------------------------------------------------
1198 ! make_ensemble1.F
1199 !--------------------------------------------------------------------------
1200       subroutine make_ensembles(islice,*)
1201 ! construct the conformational ensembles at REMD temperatures
1202       use geometry_data, only:c
1203       use io_base, only:ilen
1204       use io_wham, only:pdboutW
1205       use geometry, only:returnbox
1206 !      implicit none
1207 !      include "DIMENSIONS"
1208 !      include "DIMENSIONS.ZSCOPT"
1209 !      include "DIMENSIONS.FREE"
1210 #ifdef MPI
1211       include "mpif.h"
1212 !      include "COMMON.MPI"
1213       integer :: ierror,errcode,status(MPI_STATUS_SIZE)
1214 #endif
1215 !      include "COMMON.IOUNITS"
1216 !      include "COMMON.CONTROL"
1217 !      include "COMMON.FREE"
1218 !      include "COMMON.ENERGIES"
1219 !      include "COMMON.FFIELD"
1220 !      include "COMMON.INTERACT"
1221 !      include "COMMON.SBRIDGE"
1222 !      include "COMMON.CHAIN"
1223 !      include "COMMON.PROTFILES"
1224 !      include "COMMON.PROT"
1225       real(kind=4) :: csingle(3,nres*2)
1226       real(kind=8),dimension(6) :: fT,fTprim,fTbis
1227       real(kind=8) :: quot,quotl1,quotl,kfacl,&
1228         eprim,ebis,temper,kfac=2.4d0,T0=300.0d0
1229       real(kind=8) :: etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,&
1230             escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,&
1231             eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt, &
1232             ecation_prot, ecationcation,evdwpp,eespp,evdwpsb,eelpsb,evdwsb, &
1233             eelsb,estr_nucl,ebe_nucl,esbloc,etors_nucl,etors_d_nucl,&
1234             ecorr_nucl,ecorr3_nucl,escbase, epepbase,escpho, epeppho
1235
1236       integer :: i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
1237       real(kind=8) :: qfree,sumprob,eini,efree,rmsdev
1238       character(len=80) :: bxname
1239       character(len=2) :: licz1,licz2
1240       character(len=3) :: licz3,licz4
1241       character(len=5) :: ctemper
1242 !el      integer ilen
1243 !el      external ilen
1244       real(kind=4) :: Fdimless(MaxStr),Fdimless_(MaxStr)
1245       real(kind=8) :: enepot(MaxStr)
1246       integer :: iperm(MaxStr)
1247       integer :: islice
1248       integer,dimension(0:nprocs) :: scount_
1249
1250 #ifdef MPI
1251       if (me.eq.Master) then
1252 #endif
1253       write (licz2,'(bz,i2.2)') islice
1254       if (nslice.eq.1) then
1255         if (.not.separate_parset) then
1256           bxname = prefix(:ilen(prefix))//".bx"
1257         else
1258           write (licz3,'(bz,i3.3)') myparm
1259           bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx"
1260         endif
1261       else
1262         if (.not.separate_parset) then
1263           bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
1264         else
1265           write (licz3,'(bz,i3.3)') myparm
1266           bxname = prefix(:ilen(prefix))//"par_"//licz3// &
1267             "_slice_"//licz2//".bx"
1268         endif
1269       endif
1270       open (ientout,file=bxname,status="unknown",&
1271         form="unformatted",access="direct",recl=lenrec1)
1272 #ifdef MPI
1273       endif
1274 #endif
1275       do iparm=1,iparm
1276         if (iparm.ne.iparmprint) exit
1277         call restore_parm(iparm)
1278         do ib=1,nT_h(iparm)
1279 #ifdef DEBUG
1280           write (iout,*) "iparm",iparm," ib",ib
1281 #endif
1282           temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1283 !          quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1284 !          quotl=1.0d0
1285 !          kfacl=1.0d0
1286 !          do l=1,5
1287 !            quotl1=quotl
1288 !            quotl=quotl*quot
1289 !            kfacl=kfacl*kfac
1290 !            fT(l)=kfacl/(kfacl-1.0d0+quotl)
1291 !          enddo
1292 !el old rescale weights
1293 !
1294 !            if (rescale_mode.eq.1) then
1295 !              quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1296 #if defined(FUNCTH)
1297               tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1298               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
1299 #elif defined(FUNCT)
1300               ft(6)=quot
1301 #else
1302               ft(6)=1.0d0
1303 #endif
1304 !              quotl=1.0d0
1305 !              kfacl=1.0d0
1306 !              do l=1,5
1307 !                quotl1=quotl
1308 !                quotl=quotl*quot
1309 !                kfacl=kfacl*kfac
1310 !                fT(l)=kfacl/(kfacl-1.0d0+quotl)
1311 !              enddo
1312 !            else if (rescale_mode.eq.2) then
1313 !              quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1314 !#if defined(FUNCTH)
1315 !              tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1316 !              ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0
1317 !#elif defined(FUNCT)
1318 !              ft(6)=quot
1319 !#else
1320 !              ft(6)=1.0d0
1321 !#endif
1322 !              quotl=1.0d0
1323 !              do l=1,5
1324 !                quotl=quotl*quot
1325 !                fT(l)=1.12692801104297249644d0/ &
1326 !                   dlog(dexp(quotl)+dexp(-quotl))
1327 !              enddo
1328 !              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
1329 !            else if (rescale_mode.eq.0) then
1330 !              do l=1,5
1331 !                fT(l)=0.0d0
1332 !              enddo
1333 !            else
1334 !              write (iout,*) &
1335 !              "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
1336 !              call flush(iout)
1337 !              return 1
1338 !            endif
1339 ! el end old rescale weihgts
1340           call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3))
1341
1342 #ifdef MPI
1343           do i=1,scount(me1)
1344 #else
1345           do i=1,ntot(islice)
1346 #endif
1347             evdw=enetb(1,i,iparm)
1348 !            evdw_t=enetb(21,i,iparm)
1349             evdw_t=enetb(20,i,iparm)
1350 #ifdef SCP14
1351 !            evdw2_14=enetb(17,i,iparm)
1352             evdw2_14=enetb(18,i,iparm)
1353             evdw2=enetb(2,i,iparm)+evdw2_14
1354 #else
1355             evdw2=enetb(2,i,iparm)
1356             evdw2_14=0.0d0
1357 #endif
1358 #ifdef SPLITELE
1359             ees=enetb(3,i,iparm)
1360             evdw1=enetb(16,i,iparm)
1361 #else
1362             ees=enetb(3,i,iparm)
1363             evdw1=0.0d0
1364 #endif
1365             ecorr=enetb(4,i,iparm)
1366             ecorr5=enetb(5,i,iparm)
1367             ecorr6=enetb(6,i,iparm)
1368             eel_loc=enetb(7,i,iparm)
1369             eello_turn3=enetb(8,i,iparm)
1370             eello_turn4=enetb(9,i,iparm)
1371             eello_turn6=enetb(10,i,iparm)
1372             ebe=enetb(11,i,iparm)
1373             escloc=enetb(12,i,iparm)
1374             etors=enetb(13,i,iparm)
1375             etors_d=enetb(14,i,iparm)
1376             ehpb=enetb(15,i,iparm)
1377
1378             estr=enetb(17,i,iparm)
1379 !            estr=enetb(18,i,iparm)
1380 !            esccor=enetb(19,i,iparm)
1381             esccor=enetb(21,i,iparm)
1382 !            edihcnstr=enetb(20,i,iparm)
1383             edihcnstr=enetb(19,i,iparm)
1384             ecationcation=enetb(42,i,iparm)
1385             ecation_prot=enetb(41,i,iparm)
1386             evdwpp =      enetb(26,i,iparm)
1387             eespp =      enetb(27,i,iparm)
1388             evdwpsb =      enetb(28,i,iparm)
1389             eelpsb =      enetb(29,i,iparm)
1390             evdwsb =      enetb(30,i,iparm)
1391             eelsb =      enetb(31,i,iparm)
1392             estr_nucl =      enetb(32,i,iparm)
1393             ebe_nucl =      enetb(33,i,iparm)
1394             esbloc =      enetb(34,i,iparm)
1395             etors_nucl =      enetb(35,i,iparm)
1396             etors_d_nucl =      enetb(36,i,iparm)
1397             ecorr_nucl =      enetb(37,i,iparm)
1398             ecorr3_nucl =      enetb(38,i,iparm)
1399
1400 !#ifdef SPLITELE
1401 !            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1402 !            +wvdwpp*evdw1 &
1403 !            +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1404 !            +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1405 !            +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1406 !            +ft(2)*wturn3*eello_turn3 &
1407 !            +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc &
1408 !            +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1409 !            +wbond*estr
1410 !#else
1411 !            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1412 !            +ft(1)*welec*(ees+evdw1) &
1413 !            +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1414 !            +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1415 !            +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1416 !            +ft(2)*wturn3*eello_turn3 &
1417 !            +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1418 !            +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1419 !            +wbond*estr
1420 !#endif
1421
1422 #ifdef SPLITELE
1423             etot=wsc*evdw+wscp*evdw2+welec*ees &
1424             +wvdwpp*evdw1 &
1425             +wang*ebe+wtor*etors+wscloc*escloc &
1426             +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
1427             +wcorr6*ecorr6+wturn4*eello_turn4 &
1428             +wturn3*eello_turn3 &
1429             +wturn6*eello_turn6+wel_loc*eel_loc &
1430             +edihcnstr+wtor_d*etors_d+wsccor*esccor &
1431             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
1432             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1433             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1434             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1435             +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
1436             +wscbase*escbase&
1437             +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
1438
1439 #else
1440             etot=wsc*evdw+wscp*evdw2 &
1441             +welec*(ees+evdw1) &
1442             +wang*ebe+wtor*etors+wscloc*escloc &
1443             +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
1444             +wcorr6*ecorr6+wturn4*eello_turn4 &
1445             +wturn3*eello_turn3 &
1446             +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
1447             +wtor_d*etors_d+wsccor*esccor &
1448             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
1449             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1450             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1451             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1452             +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
1453             +wscbase*escbase&
1454             +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
1455
1456 #endif
1457
1458 #ifdef MPI
1459             Fdimless(i)= &
1460               beta_h(ib,iparm)*etot-entfac(i)
1461             potE(i,iparm)=etot
1462 #ifdef DEBUG
1463             write (iout,*) i,indstart(me)+i-1,ib,&
1464              1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),&
1465              -entfac(i),Fdimless(i)
1466 #endif
1467 #else
1468             Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
1469             potE(i,iparm)=etot
1470 #endif
1471           enddo   ! i
1472 #ifdef MPI
1473           do i=1,scount(me1)
1474             Fdimless_(i)=Fdimless(i)
1475           enddo
1476           call MPI_Gatherv(Fdimless_(1),scount(me),&
1477            MPI_REAL,Fdimless(1),&
1478            scount_(0),idispl(0),MPI_REAL,Master,&
1479            WHAM_COMM, IERROR)
1480 #ifdef DEBUG
1481           call MPI_Gatherv(potE(1,iparm),scount_(me),&
1482            MPI_DOUBLE_PRECISION,potE(1,iparm),&
1483            scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1484            WHAM_COMM, IERROR)
1485           call MPI_Gatherv(entfac(1),scount(me),&
1486            MPI_DOUBLE_PRECISION,entfac(1),&
1487            scount_(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1488            WHAM_COMM, IERROR)
1489 #endif
1490           if (me.eq.Master) then
1491 #ifdef DEBUG
1492           write (iout,*) "The FDIMLESS array before sorting"
1493           do i=1,ntot(islice)
1494             write (iout,*) i,fdimless(i)
1495           enddo
1496 #endif
1497 #endif
1498           do i=1,ntot(islice)
1499             iperm(i)=i
1500           enddo
1501           call mysort1(ntot(islice),Fdimless,iperm)
1502 #ifdef DEBUG
1503           write (iout,*) "The FDIMLESS array after sorting"
1504           do i=1,ntot(islice)
1505             write (iout,*) i,iperm(i),fdimless(i)
1506           enddo
1507 #endif
1508           qfree=0.0d0
1509           do i=1,ntot(islice)
1510             qfree=qfree+exp(-fdimless(i)+fdimless(1))
1511           enddo
1512 !          write (iout,*) "qfree",qfree
1513           nlist=1
1514           sumprob=0.0
1515           do i=1,min0(ntot(islice),ensembles)
1516             sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
1517 #ifdef DEBUG
1518             write (iout,*) i,ib,beta_h(ib,iparm),&
1519              1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),&
1520              potE(iperm(i),iparm),&
1521              -entfac(iperm(i)),fdimless(i),sumprob
1522 #endif
1523             if (sumprob.gt.0.99d0) goto 122
1524             nlist=nlist+1
1525           enddo
1526   122     continue
1527 #ifdef MPI
1528           endif
1529           call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,&
1530              IERROR)
1531           call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,&
1532              IERROR)
1533           do i=1,nlist
1534             ii=iperm(i)
1535             iproc=0
1536             do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
1537               iproc=iproc+1
1538             enddo
1539             if (iproc.ge.nprocs) then
1540               write (iout,*) "Fatal error: processor out of range",iproc
1541               call flush(iout)
1542               if (bxfile) then
1543                 close (ientout)
1544               else
1545                 close (ientout,status="delete")
1546               endif
1547               return 1
1548             endif
1549             ik=ii-indstart(iproc)+1
1550             if (iproc.ne.Master) then
1551               if (me.eq.iproc) then
1552 #ifdef DEBUG
1553                 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,&
1554                  " energy",potE(ik,iparm)
1555 #endif
1556                 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,&
1557                   Master,i,WHAM_COMM,IERROR)
1558               else if (me.eq.Master) then
1559                 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,&
1560                   WHAM_COMM,STATUS,IERROR)
1561               endif
1562             else if (me.eq.Master) then
1563               enepot(i)=potE(ik,iparm)
1564             endif
1565           enddo
1566 #else
1567           do i=1,nlist
1568             enepot(i)=potE(iperm(i),iparm)
1569           enddo
1570 #endif
1571 #ifdef MPI
1572           if (me.eq.Master) then
1573 #endif
1574           write(licz3,'(bz,i3.3)') iparm
1575           write(licz2,'(bz,i2.2)') islice
1576           if (temper.lt.100.0d0) then
1577             write(ctemper,'(f3.0)') temper
1578           else if (temper.lt.1000.0) then
1579             write (ctemper,'(f4.0)') temper
1580           else
1581             write (ctemper,'(f5.0)') temper
1582           endif
1583           if (nparmset.eq.1) then
1584             if (separate_parset) then
1585               write(licz4,'(bz,i3.3)') myparm
1586               pdbname=prefix(:ilen(prefix))//"_par"//licz4
1587             else
1588               pdbname=prefix(:ilen(prefix))
1589             endif
1590           else
1591             pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
1592           endif
1593           if (nslice.eq.1) then
1594             pdbname=pdbname(:ilen(pdbname))//"_T_"// &
1595               ctemper(:ilen(ctemper))//"pdb"
1596           else
1597             pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"// &
1598               ctemper(:ilen(ctemper))//"pdb"
1599           endif
1600           open(ipdb,file=pdbname)
1601           do i=1,nlist
1602             read (ientout,rec=iperm(i)) &
1603               ((csingle(l,k),l=1,3),k=1,nres),&
1604               ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1605               nss,(ihpb(k),jhpb(k),k=1,nss),&
1606               eini,efree,rmsdev,iscor
1607             do j=1,2*nres
1608               do k=1,3
1609                 c(k,j)=csingle(k,j)
1610               enddo
1611             enddo
1612             call returnbox
1613             do j=1,2*nres
1614               do k=1,3
1615                 csingle(k,j)=c(k,j)
1616               enddo
1617             enddo
1618
1619             eini=fdimless(i)
1620             call pdboutW(iperm(i),temper,eini,enepot(i),efree,rmsdev)
1621           enddo
1622 #ifdef MPI
1623           endif
1624 #endif
1625         enddo     ! ib
1626       enddo       ! iparm
1627       if (bxfile) then
1628         close(ientout)
1629       else
1630         close(ientout,status="delete")
1631       endif
1632       do i=0,nprocs
1633         scount(i)=scount_(i)
1634       enddo
1635       return
1636       end subroutine make_ensembles
1637 !--------------------------------------------------------------------------
1638       subroutine mysort1(n, x, ipermut)
1639 !      implicit none
1640       integer :: i,j,imax,ipm,n
1641       real(kind=4) :: x(n)
1642       integer :: ipermut(n)
1643       real(kind=4) :: xtemp
1644       do i=1,n
1645         xtemp=x(i)
1646         imax=i
1647         do j=i+1,n
1648           if (x(j).lt.xtemp) then
1649             imax=j
1650             xtemp=x(j)
1651           endif
1652         enddo
1653         x(imax)=x(i)
1654         x(i)=xtemp
1655         ipm=ipermut(imax)
1656         ipermut(imax)=ipermut(i)
1657         ipermut(i)=ipm
1658       enddo
1659       return
1660       end subroutine mysort1
1661 !--------------------------------------------------------------------------
1662       subroutine alloc_enecalc_arrays(iparm)
1663
1664       use control_data
1665       use geometry_data, only:maxlob
1666       integer :: iparm
1667 !---------------------------
1668 ! COMMON.ENERGIES form wham_data
1669 !      common /energies/
1670       allocate(potE(MaxStr_Proc,iparm)) !(MaxStr_Proc,Max_Parm)
1671       allocate(entfac(MaxStr_Proc)) !(MaxStr_Proc)
1672       allocate(q(nQ+2,MaxStr_Proc)) !(MaxQ+2,MaxStr_Proc)
1673       allocate(enetb(max_ene,MaxStr_Proc,iparm)) !(max_ene,MaxStr_Proc,Max_Parm)
1674 !
1675 ! allocate ENECALC arrays
1676 !---------------------------
1677 ! COMMON.ALLPARM
1678 !      common /allparm/
1679       allocate(ww_all(max_eneW,iparm)) !(max_ene,max_parm) ! max_eneW
1680       allocate(vbldp0_all(iparm),akp_all(nParmSet)) !(max_parm)
1681       allocate(vbldsc0_all(maxbondterm,ntyp,iparm),&
1682         aksc_all(maxbondterm,ntyp,iparm),&
1683         abond0_all(maxbondterm,ntyp,iparm)) !(maxbondterm,ntyp,max_parm)
1684       allocate(a0thet_all(-ntyp:ntyp,iparm)) !(-ntyp:ntyp,max_parm)
1685       allocate(athet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm),&
1686         bthet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm)) !(2,-ntyp:ntyp,-1:1,-1:1,max_parm)
1687       allocate(polthet_all(0:3,-ntyp:ntyp,iparm)) !(0:3,-ntyp:ntyp,max_parm)
1688       allocate(gthet_all(3,-ntyp:ntyp,iparm)) !(3,-ntyp:ntyp,max_parm)
1689       allocate(theta0_all(-ntyp:ntyp,iparm),&
1690         sig0_all(-ntyp:ntyp,iparm),sigc0_all(-ntyp:ntyp,nParmSet)) !(-ntyp:ntyp,max_parm)
1691       allocate(aa0thet_all(-nthetyp-1:nthetyp+1,&
1692         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1693 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1694       allocate(aathet_all(ntheterm,-nthetyp-1:nthetyp+1,&
1695         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1696 !(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1697       allocate(bbthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1698         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1699       allocate(ccthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1700         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1701       allocate(ddthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1702         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1703       allocate(eethet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1704         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1705 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,
1706 !     & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1707       allocate(ffthet_all1(ndouble,ndouble,ntheterm3,&
1708         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1709         -nthetyp-1:nthetyp+1,iparm))
1710       allocate(ggthet_all1(ndouble,ndouble,ntheterm3,&
1711         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1712         -nthetyp-1:nthetyp+1,iparm))
1713       allocate(ffthet_all2(ndouble,ndouble,ntheterm3,&
1714         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1715         -nthetyp-1:nthetyp+1,iparm))
1716       allocate(ggthet_all2(ndouble,ndouble,ntheterm3,&
1717         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1718         -nthetyp-1:nthetyp+1,iparm))
1719 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1720 !-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm)
1721       allocate(dsc_all(ntyp1,iparm),dsc0_all(ntyp1,nParmSet)) !(ntyp1,max_parm)
1722       allocate(bsc_all(maxlob,ntyp,iparm))
1723 !(maxlob,ntyp,max_parm)
1724       allocate(censc_all(3,maxlob,-ntyp:ntyp,iparm)) !(3,maxlob,-ntyp:ntyp,max_parm)
1725       allocate(gaussc_all(3,3,maxlob,-ntyp:ntyp,iparm)) !(3,3,maxlob,-ntyp:ntyp,max_parm)
1726       allocate(sc_parmin_all(65,ntyp,iparm)) !(65,ntyp,max_parm)
1727       allocate(v0_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1728 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1729       allocate(v1_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1730       allocate(v2_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1731 !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1732       allocate(vlor1_all(maxlor,ntortyp,ntortyp,iparm))
1733       allocate(vlor2_all(maxlor,ntortyp,ntortyp,iparm))
1734       allocate(vlor3_all(maxlor,ntortyp,ntortyp,iparm)) !(maxlor,maxtor,maxtor,max_parm)
1735       allocate(v1c_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1736         -ntortyp:ntortyp,2,iparm))
1737       allocate(v1s_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1738         -ntortyp:ntortyp,2,iparm))
1739 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1740       allocate(v2c_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1741         -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1742 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1743       allocate(v2s_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1744         -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1745 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1746       allocate(b1_all(2,-nloctyp:nloctyp,iparm))
1747       allocate(b2_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1748       allocate(cc_all(2,2,-nloctyp:nloctyp,iparm))
1749       allocate(dd_all(2,2,-nloctyp:nloctyp,iparm))
1750       allocate(ee_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1751       allocate(ctilde_all(2,2,-nloctyp:nloctyp,iparm))
1752       allocate(dtilde_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1753       allocate(b1tilde_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1754       allocate(app_all(2,2,iparm),bpp_all(2,2,nParmSet),&
1755         ael6_all(2,2,iparm),ael3_all(2,2,nParmSet)) !(2,2,max_parm)
1756       allocate(aad_all(ntyp,2,iparm),bad_all(ntyp,2,nParmSet)) !(ntyp,2,max_parm)
1757       allocate(aa_all(ntyp,ntyp,iparm),bb_all(ntyp,ntyp,nParmSet),&
1758         augm_all(ntyp,ntyp,iparm),eps_all(ntyp,ntyp,nParmSet),&
1759         sigma_all(ntyp,ntyp,iparm),r0_all(ntyp,ntyp,nParmSet),&
1760         chi_all(ntyp,ntyp,iparm)) !(ntyp,ntyp,max_parm)
1761       allocate(chip_all(ntyp,iparm),alp_all(ntyp,nParmSet)) !(ntyp,max_parm)
1762       allocate(ebr_all(iparm),d0cm_all(nParmSet),akcm_all(nParmSet),&
1763         akth_all(iparm),akct_all(nParmSet),v1ss_all(nParmSet),&
1764         v2ss_all(iparm),v3ss_all(nParmSet)) !(max_parm)
1765       allocate(v1sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1766       allocate(v2sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1767 !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
1768       allocate(nlob_all(ntyp1,iparm)) !(ntyp1,max_parm)
1769       allocate(nlor_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1770       allocate(nterm_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1771 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1772       allocate(ntermd1_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1773         -ntortyp:ntortyp,2,iparm))
1774       allocate(ntermd2_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1775         -ntortyp:ntortyp,2,iparm))
1776 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1777       allocate(nbondterm_all(ntyp,iparm)) !(ntyp,max_parm)
1778       allocate(ithetyp_all(-ntyp1:ntyp1,iparm)) !(-ntyp1:ntyp1,max_parm)
1779       allocate(nthetyp_all(iparm),ntheterm_all(nParmSet),&
1780         ntheterm2_all(iparm),ntheterm3_all(nParmSet),&
1781         nsingle_all(iparm),&
1782         ndouble_all(iparm),nntheterm_all(nParmSet)) !(max_parm)
1783       allocate(nterm_sccor_all(-ntyp:ntyp,-ntyp:ntyp,iparm)) !(-ntyp:ntyp,-ntyp:ntyp,max_parm)
1784 !
1785       end subroutine alloc_enecalc_arrays
1786 !--------------------------------------------------------------------------
1787 !--------------------------------------------------------------------------
1788       end module ene_calc