source/wham/enecalc.f90

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