source/wham/enecalc.F90

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