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,49)
 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,49
 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,mnum,mnum1
 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         mnum=molnum(j)
 541         mnum1=molnum(j-1)
 542         if (mnum.ne.1) cycle
 543         if (itype(j-1,mnum1).ne.ntyp1_molec(mnum1) .and. itype(j,mnum).ne.ntyp1_molec(mnum) .and. &
 544           (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
 545           if (iprint.gt.0) &
 546           write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),&
 547             " for conformation",ii
 548           if (iprint.gt.1) then
 549             write (iout,*) "The Cartesian geometry is:"
 550             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 551             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 552             write (iout,*) "The internal geometry is:"
 553             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 554             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 555             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 556             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 557             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 558             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 559           endif
 560           if (iprint.gt.0) write (iout,*) &
 561             "This conformation WILL NOT be added to the database."
 562           conf_check=.false.
 563           return
 564         endif
 565       enddo
 566       do j=nnt,nct
 567         mnum=molnum(j)
 568        if (mnum.gt.3) cycle
 569         itj=itype(j,mnum)
 570         if (itype(j,mnum).ne.10 .and.itype(j,mnum).ne.ntyp1 .and. &
 571            (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then
 572           if (iprint.gt.0) &
 573           write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),&
 574            " for conformation",ii
 575           if (iprint.gt.1) then
 576             write (iout,*) "The Cartesian geometry is:"
 577             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 578             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 579             write (iout,*) "The internal geometry is:"
 580             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 581             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 582             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 583             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 584             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 585             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 586           endif
 587           if (iprint.gt.0) write (iout,*) &
 588             "This conformation WILL NOT be added to the database."
 589           conf_check=.false.
 590           return
 591         endif
 592       enddo
 593       do j=3,nres
 594         if (theta(j).le.0.0d0) then
 595           if (iprint.gt.0) &
 596           write (iout,*) "Zero theta angle(s) in conformation",ii
 597           if (iprint.gt.1) then
 598             write (iout,*) "The Cartesian geometry is:"
 599             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 600             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 601             write (iout,*) "The internal geometry is:"
 602             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 603             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 604             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 605             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 606             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 607             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 608           endif
 609           if (iprint.gt.0) write (iout,*) &
 610             "This conformation WILL NOT be added to the database."
 611           conf_check=.false.
 612           return
 613         endif
 614         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
 615       enddo
 616       conf_check=.true.
 617 !      write (iout,*) "conf_check passed",ii
 618       return
 619       end function conf_check
 620 !-----------------------------------------------------------------------------
 621 ! store_parm.F
 622 !-----------------------------------------------------------------------------
 623       subroutine store_parm(iparm)
 624 !
 625 ! Store parameters of set IPARM
 626 ! valence angles and the side chains and energy parameters.
 627 !
 628 !      implicit none
 629 !      include 'DIMENSIONS'
 630 !      include 'DIMENSIONS.ZSCOPT'
 631 !      include 'DIMENSIONS.FREE'
 632 !      include 'COMMON.IOUNITS'
 633 !      include 'COMMON.CHAIN'
 634 !      include 'COMMON.INTERACT'
 635 !      include 'COMMON.GEO'
 636 !      include 'COMMON.LOCAL'
 637 !      include 'COMMON.TORSION'
 638 !      include 'COMMON.FFIELD'
 639 !      include 'COMMON.NAMES'
 640 !      include 'COMMON.SBRIDGE'
 641 !      include 'COMMON.SCROT'
 642 !      include 'COMMON.SCCOR'
 643 !      include 'COMMON.ALLPARM'
 644       integer :: i,j,k,l,m,mm,iparm,ichir1,ichir2,iblock,iii
 645
 646       call alloc_enecalc_arrays(iparm)
 647 !el      allocate(ww_all(n_ene,iparm))
 648 ! Store weights
 649       ww_all(1,iparm)=wsc
 650       ww_all(2,iparm)=wscp
 651       ww_all(3,iparm)=welec
 652       ww_all(4,iparm)=wcorr
 653       ww_all(5,iparm)=wcorr5
 654       ww_all(6,iparm)=wcorr6
 655       ww_all(7,iparm)=wel_loc
 656       ww_all(8,iparm)=wturn3
 657       ww_all(9,iparm)=wturn4
 658       ww_all(10,iparm)=wturn6
 659       ww_all(11,iparm)=wang
 660       ww_all(12,iparm)=wscloc
 661       ww_all(13,iparm)=wtor
 662       ww_all(14,iparm)=wtor_d
 663       ww_all(15,iparm)=wstrain
 664       ww_all(16,iparm)=wvdwpp
 665       ww_all(17,iparm)=wbond
 666       ww_all(19,iparm)=wsccor
 667 ! Store bond parameters
 668       vbldp0_all(iparm)=vbldp0
 669       akp_all(iparm)=akp
 670       do i=1,ntyp
 671         nbondterm_all(i,iparm)=nbondterm(i)
 672         do j=1,nbondterm(i)
 673           vbldsc0_all(j,i,iparm)=vbldsc0(j,i)
 674           aksc_all(j,i,iparm)=aksc(j,i)
 675           abond0_all(j,i,iparm)=abond0(j,i)
 676         enddo
 677       enddo
 678 ! Store bond angle parameters
 679 #ifdef CRYST_THETA
 680       do i=-ntyp,ntyp
 681         a0thet_all(i,iparm)=a0thet(i)
 682         do ichir1=-1,1
 683         do ichir2=-1,1
 684         do j=1,2
 685           athet_all(j,i,ichir1,ichir2,iparm)=athet(j,i,ichir1,ichir2)
 686           bthet_all(j,i,ichir1,ichir2,iparm)=bthet(j,i,ichir1,ichir2)
 687         enddo
 688         enddo
 689         enddo
 690         do j=0,3
 691           polthet_all(j,i,iparm)=polthet(j,i)
 692         enddo
 693         do j=1,3
 694           gthet_all(j,i,iparm)=gthet(j,i)
 695         enddo
 696         theta0_all(i,iparm)=theta0(i)
 697         sig0_all(i,iparm)=sig0(i)
 698         sigc0_all(i,iparm)=sigc0(i)
 699       enddo
 700 #else
 701       nthetyp_all(iparm)=nthetyp
 702       ntheterm_all(iparm)=ntheterm
 703       ntheterm2_all(iparm)=ntheterm2
 704       ntheterm3_all(iparm)=ntheterm3
 705       nsingle_all(iparm)=nsingle
 706       ndouble_all(iparm)=ndouble
 707       nntheterm_all(iparm)=nntheterm
 708       do i=-ntyp,ntyp
 709         ithetyp_all(i,iparm)=ithetyp(i)
 710       enddo
 711       do iblock=1,2
 712       do i=-nthetyp-1,nthetyp+1
 713         do j=-nthetyp-1,nthetyp+1
 714           do k=-nthetyp-1,nthetyp+1
 715             aa0thet_all(i,j,k,iblock,iparm)=aa0thet(i,j,k,iblock)
 716             do l=1,ntheterm
 717               aathet_all(l,i,j,k,iblock,iparm)=aathet(l,i,j,k,iblock)
 718             enddo
 719             do l=1,ntheterm2
 720               do m=1,nsingle
 721                 bbthet_all(m,l,i,j,k,iblock,iparm)= &
 722       bbthet(m,l,i,j,k,iblock)
 723                 ccthet_all(m,l,i,j,k,iblock,iparm)= &
 724       ccthet(m,l,i,j,k,iblock)
 725                 ddthet_all(m,l,i,j,k,iblock,iparm)= &
 726       ddthet(m,l,i,j,k,iblock)
 727                 eethet_all(m,l,i,j,k,iblock,iparm)= &
 728       eethet(m,l,i,j,k,iblock)
 729               enddo
 730             enddo
 731             do l=1,ntheterm3
 732               do m=1,ndouble
 733                 do mm=1,ndouble
 734                 if (iblock.eq.1) then
 735                  ffthet_all1(mm,m,l,i,j,k,iparm)=&
 736          ffthet(mm,m,l,i,j,k,iblock)
 737                  ggthet_all1(mm,m,l,i,j,k,iparm)=&
 738       ggthet(mm,m,l,i,j,k,iblock)
 739                   else
 740                  ffthet_all2(mm,m,l,i,j,k,iparm)=&
 741          ffthet(mm,m,l,i,j,k,iblock)
 742                  ggthet_all2(mm,m,l,i,j,k,iparm)=&
 743       ggthet(mm,m,l,i,j,k,iblock)
 744                   endif
 745                 enddo
 746               enddo
 747             enddo
 748           enddo
 749         enddo
 750       enddo
 751       enddo
 752 #endif
 753 #ifdef CRYST_SC
 754 ! Store the sidechain rotamer parameters
 755       do i=-ntyp,ntyp
 756        iii=iabs(i)
 757 !!       write (iout,*) i,"storeparm1"
 758        if (i.eq.0) cycle
 759         nlob_all(iii,iparm)=nlob(iii)
 760         do j=1,nlob(iii)
 761           bsc_all(j,iii,iparm)=bsc(j,iii)
 762           do k=1,3
 763             censc_all(k,j,i,iparm)=censc(k,j,i)
 764           enddo
 765           do k=1,3
 766             do l=1,3
 767               gaussc_all(l,k,j,i,iparm)=gaussc(l,k,j,i)
 768             enddo
 769           enddo
 770         enddo
 771       enddo
 772 #else
 773       do i=1,ntyp
 774         do j=1,65
 775           sc_parmin_all(j,i,iparm)=sc_parmin(j,i)
 776         enddo
 777       enddo
 778 #endif
 779 ! Store the torsional parameters
 780       do iblock=1,2
 781       do i=-ntortyp+1,ntortyp-1
 782         do j=-ntortyp+1,ntortyp-1
 783           v0_all(i,j,iblock,iparm)=v0(i,j,iblock)
 784           nterm_all(i,j,iblock,iparm)=nterm(i,j,iblock)
 785           nlor_all(i,j,iblock,iparm)=nlor(i,j,iblock)
 786           do k=1,nterm(i,j,iblock)
 787             v1_all(k,i,j,iblock,iparm)=v1(k,i,j,iblock)
 788             v2_all(k,i,j,iblock,iparm)=v2(k,i,j,iblock)
 789           enddo
 790           do k=1,nlor(i,j,iblock)
 791             vlor1_all(k,i,j,iparm)=vlor1(k,i,j)
 792             vlor2_all(k,i,j,iparm)=vlor2(k,i,j)
 793             vlor3_all(k,i,j,iparm)=vlor3(k,i,j)
 794           enddo
 795         enddo
 796       enddo
 797       enddo
 798 ! Store the double torsional parameters
 799       do iblock=1,2
 800       do i=-ntortyp+1,ntortyp-1
 801         do j=-ntortyp+1,ntortyp-1
 802           do k=-ntortyp+1,ntortyp-1
 803             ntermd1_all(i,j,k,iblock,iparm)=ntermd_1(i,j,k,iblock)
 804             ntermd2_all(i,j,k,iblock,iparm)=ntermd_2(i,j,k,iblock)
 805             do l=1,ntermd_1(i,j,k,iblock)
 806               v1c_all(1,l,i,j,k,iblock,iparm)=v1c(1,l,i,j,k,iblock)
 807               v1c_all(2,l,i,j,k,iblock,iparm)=v1c(2,l,i,j,k,iblock)
 808               v2c_all(1,l,i,j,k,iblock,iparm)=v2c(1,l,i,j,k,iblock)
 809               v2c_all(2,l,i,j,k,iblock,iparm)=v2c(2,l,i,j,k,iblock)
 810             enddo
 811             do l=1,ntermd_2(i,j,k,iblock)
 812               do m=1,ntermd_2(i,j,k,iblock)
 813                 v2s_all(l,m,i,j,k,iblock,iparm)=v2s(l,m,i,j,k,iblock)
 814               enddo
 815             enddo
 816           enddo
 817         enddo
 818       enddo
 819       enddo
 820 ! Store parameters of the cumulants
 821       do i=-nloctyp,nloctyp
 822         do j=1,2
 823           b1_all(j,i,iparm)=b1(j,i)
 824           b1tilde_all(j,i,iparm)=b1tilde(j,i)
 825           b2_all(j,i,iparm)=b2(j,i)
 826         enddo
 827         do j=1,2
 828           do k=1,2
 829             cc_all(k,j,i,iparm)=cc(k,j,i)
 830             ctilde_all(k,j,i,iparm)=ctilde(k,j,i)
 831             dd_all(k,j,i,iparm)=dd(k,j,i)
 832             dtilde_all(k,j,i,iparm)=dtilde(k,j,i)
 833             ee_all(k,j,i,iparm)=ee(k,j,i)
 834           enddo
 835         enddo
 836       enddo
 837 ! Store the parameters of electrostatic interactions
 838       do i=1,2
 839         do j=1,2
 840           app_all(j,i,iparm)=app(j,i)
 841           bpp_all(j,i,iparm)=bpp(j,i)
 842           ael6_all(j,i,iparm)=ael6(j,i)
 843           ael3_all(j,i,iparm)=ael3(j,i)
 844         enddo
 845       enddo
 846 ! Store sidechain parameters
 847       do i=1,ntyp
 848         do j=1,ntyp
 849           aa_all(j,i,iparm)=aa_aq(j,i)
 850           bb_all(j,i,iparm)=bb_aq(j,i)
 851           r0_all(j,i,iparm)=r0(j,i)
 852           sigma_all(j,i,iparm)=sigma(j,i)
 853           chi_all(j,i,iparm)=chi(j,i)
 854           augm_all(j,i,iparm)=augm(j,i)
 855           eps_all(j,i,iparm)=eps(j,i)
 856         enddo
 857       enddo
 858       do i=1,ntyp
 859         chip_all(i,iparm)=chip(i)
 860         alp_all(i,iparm)=alp(i)
 861       enddo
 862 ! Store the SCp parameters
 863       do i=1,ntyp
 864         do j=1,2
 865           aad_all(i,j,iparm)=aad(i,j)
 866           bad_all(i,j,iparm)=bad(i,j)
 867         enddo
 868       enddo
 869 ! Store disulfide-bond parameters
 870       ebr_all(iparm)=ebr
 871       d0cm_all(iparm)=d0cm
 872       akcm_all(iparm)=akcm
 873       akth_all(iparm)=akth
 874       akct_all(iparm)=akct
 875       v1ss_all(iparm)=v1ss
 876       v2ss_all(iparm)=v2ss
 877       v3ss_all(iparm)=v3ss
 878 ! Store SC-backbone correlation parameters
 879       do i=-nsccortyp,nsccortyp
 880        do j=-nsccortyp,nsccortyp
 881
 882       nterm_sccor_all(j,i,iparm)=nterm_sccor(j,i)
 883 !      do i=1,20
 884 !        do j=1,20
 885          do l=1,3
 886           do k=1,nterm_sccor(j,i)
 887             v1sccor_all(k,l,j,i,iparm)=v1sccor(k,l,j,i)
 888             v2sccor_all(k,l,j,i,iparm)=v2sccor(k,l,j,i)
 889            enddo
 890           enddo
 891         enddo
 892       enddo
 893 write(iout,*)"end of store_parm"
 894       return
 895       end subroutine store_parm
 896 !--------------------------------------------------------------------------
 897       subroutine restore_parm(iparm)
 898 !
 899 ! Store parameters of set IPARM
 900 ! valence angles and the side chains and energy parameters.
 901 !
 902 !      implicit none
 903 !      include 'DIMENSIONS'
 904 !      include 'DIMENSIONS.ZSCOPT'
 905 !      include 'DIMENSIONS.FREE'
 906 !      include 'COMMON.IOUNITS'
 907 !      include 'COMMON.CHAIN'
 908 !      include 'COMMON.INTERACT'
 909 !      include 'COMMON.GEO'
 910 !      include 'COMMON.LOCAL'
 911 !      include 'COMMON.TORSION'
 912 !      include 'COMMON.FFIELD'
 913 !      include 'COMMON.NAMES'
 914 !      include 'COMMON.SBRIDGE'
 915 !      include 'COMMON.SCROT'
 916 !      include 'COMMON.SCCOR'
 917 !      include 'COMMON.ALLPARM'
 918       integer :: i,j,k,l,m,mm,iparm,ichir1,ichir2,iblock,iii
 919
 920 ! Restore weights
 921       wsc=ww_all(1,iparm)
 922       wscp=ww_all(2,iparm)
 923       welec=ww_all(3,iparm)
 924       wcorr=ww_all(4,iparm)
 925       wcorr5=ww_all(5,iparm)
 926       wcorr6=ww_all(6,iparm)
 927       wel_loc=ww_all(7,iparm)
 928       wturn3=ww_all(8,iparm)
 929       wturn4=ww_all(9,iparm)
 930       wturn6=ww_all(10,iparm)
 931       wang=ww_all(11,iparm)
 932       wscloc=ww_all(12,iparm)
 933       wtor=ww_all(13,iparm)
 934       wtor_d=ww_all(14,iparm)
 935       wstrain=ww_all(15,iparm)
 936       wvdwpp=ww_all(16,iparm)
 937       wbond=ww_all(17,iparm)
 938       wsccor=ww_all(19,iparm)
 939 ! Restore bond parameters
 940       vbldp0=vbldp0_all(iparm)
 941       akp=akp_all(iparm)
 942       do i=1,ntyp
 943         nbondterm(i)=nbondterm_all(i,iparm)
 944         do j=1,nbondterm(i)
 945           vbldsc0(j,i)=vbldsc0_all(j,i,iparm)
 946           aksc(j,i)=aksc_all(j,i,iparm)
 947           abond0(j,i)=abond0_all(j,i,iparm)
 948         enddo
 949       enddo
 950 ! Restore bond angle parameters
 951 #ifdef CRYST_THETA
 952       do i=-ntyp,ntyp
 953         a0thet(i)=a0thet_all(i,iparm)
 954         do ichir1=-1,1
 955         do ichir2=-1,1
 956         do j=1,2
 957           athet(j,i,ichir1,ichir2)=athet_all(j,i,ichir1,ichir2,iparm)
 958           bthet(j,i,ichir1,ichir2)=bthet_all(j,i,ichir1,ichir2,iparm)
 959         enddo
 960         enddo
 961         enddo
 962         do j=0,3
 963           polthet(j,i)=polthet_all(j,i,iparm)
 964         enddo
 965         do j=1,3
 966           gthet(j,i)=gthet_all(j,i,iparm)
 967         enddo
 968         theta0(i)=theta0_all(i,iparm)
 969         sig0(i)=sig0_all(i,iparm)
 970         sigc0(i)=sigc0_all(i,iparm)
 971       enddo
 972 #else
 973       nthetyp=nthetyp_all(iparm)
 974       ntheterm=ntheterm_all(iparm)
 975       ntheterm2=ntheterm2_all(iparm)
 976       ntheterm3=ntheterm3_all(iparm)
 977       nsingle=nsingle_all(iparm)
 978       ndouble=ndouble_all(iparm)
 979       nntheterm=nntheterm_all(iparm)
 980       do i=-ntyp,ntyp
 981         ithetyp(i)=ithetyp_all(i,iparm)
 982       enddo
 983       do iblock=1,2
 984       do i=-nthetyp-1,nthetyp+1
 985         do j=-nthetyp-1,nthetyp+1
 986           do k=-nthetyp-1,nthetyp+1
 987             aa0thet(i,j,k,iblock)=aa0thet_all(i,j,k,iblock,iparm)
 988             do l=1,ntheterm
 989               aathet(l,i,j,k,iblock)=aathet_all(l,i,j,k,iblock,iparm)
 990             enddo
 991             do l=1,ntheterm2
 992               do m=1,nsingle
 993                 bbthet(m,l,i,j,k,iblock)= &
 994       bbthet_all(m,l,i,j,k,iblock,iparm)
 995                 ccthet(m,l,i,j,k,iblock)= &
 996       ccthet_all(m,l,i,j,k,iblock,iparm)
 997                 ddthet(m,l,i,j,k,iblock)= &
 998       ddthet_all(m,l,i,j,k,iblock,iparm)
 999                 eethet(m,l,i,j,k,iblock)= &
1000       eethet_all(m,l,i,j,k,iblock,iparm)
1001               enddo
1002             enddo
1003             do l=1,ntheterm3
1004               do m=1,ndouble
1005                 do mm=1,ndouble
1006                 if (iblock.eq.1) then
1007                  ffthet(mm,m,l,i,j,k,iblock)= &
1008       ffthet_all1(mm,m,l,i,j,k,iparm)
1009                  ggthet(mm,m,l,i,j,k,iblock)= &
1010       ggthet_all1(mm,m,l,i,j,k,iparm)
1011                 else
1012                  ffthet(mm,m,l,i,j,k,iblock)= &
1013       ffthet_all2(mm,m,l,i,j,k,iparm)
1014                  ggthet(mm,m,l,i,j,k,iblock)= &
1015       ggthet_all2(mm,m,l,i,j,k,iparm)
1016                 endif
1017                 enddo
1018               enddo
1019             enddo
1020           enddo
1021         enddo
1022       enddo
1023       enddo
1024 #endif
1025 ! Restore the sidechain rotamer parameters
1026 #ifdef CRYST_SC
1027       do i=-ntyp,ntyp
1028         if (i.eq.0) cycle
1029         iii=iabs(i)
1030         nlob(iii)=nlob_all(iii,iparm)
1031         do j=1,nlob(iii)
1032           bsc(j,iii)=bsc_all(j,iii,iparm)
1033           do k=1,3
1034             censc(k,j,i)=censc_all(k,j,i,iparm)
1035           enddo
1036           do k=1,3
1037             do l=1,3
1038               gaussc(l,k,j,i)=gaussc_all(l,k,j,i,iparm)
1039             enddo
1040           enddo
1041         enddo
1042       enddo
1043 #else
1044       do i=1,ntyp
1045         do j=1,65
1046           sc_parmin(j,i)=sc_parmin_all(j,i,iparm)
1047         enddo
1048       enddo
1049 #endif
1050 ! Restore the torsional parameters
1051       do iblock=1,2
1052       do i=-ntortyp+1,ntortyp-1
1053         do j=-ntortyp+1,ntortyp-1
1054           v0(i,j,iblock)=v0_all(i,j,iblock,iparm)
1055           nterm(i,j,iblock)=nterm_all(i,j,iblock,iparm)
1056           nlor(i,j,iblock)=nlor_all(i,j,iblock,iparm)
1057           do k=1,nterm(i,j,iblock)
1058             v1(k,i,j,iblock)=v1_all(k,i,j,iblock,iparm)
1059             v2(k,i,j,iblock)=v2_all(k,i,j,iblock,iparm)
1060           enddo
1061           do k=1,nlor(i,j,iblock)
1062             vlor1(k,i,j)=vlor1_all(k,i,j,iparm)
1063             vlor2(k,i,j)=vlor2_all(k,i,j,iparm)
1064             vlor3(k,i,j)=vlor3_all(k,i,j,iparm)
1065           enddo
1066         enddo
1067       enddo
1068       enddo
1069 ! Restore the double torsional parameters
1070       do iblock=1,2
1071       do i=-ntortyp+1,ntortyp-1
1072         do j=-ntortyp+1,ntortyp-1
1073           do k=-ntortyp+1,ntortyp-1
1074             ntermd_1(i,j,k,iblock)=ntermd1_all(i,j,k,iblock,iparm)
1075             ntermd_2(i,j,k,iblock)=ntermd2_all(i,j,k,iblock,iparm)
1076             do l=1,ntermd_1(i,j,k,iblock)
1077               v1c(1,l,i,j,k,iblock)=v1c_all(1,l,i,j,k,iblock,iparm)
1078               v1c(2,l,i,j,k,iblock)=v1c_all(2,l,i,j,k,iblock,iparm)
1079               v2c(1,l,i,j,k,iblock)=v2c_all(1,l,i,j,k,iblock,iparm)
1080               v2c(2,l,i,j,k,iblock)=v2c_all(2,l,i,j,k,iblock,iparm)
1081             enddo
1082             do l=1,ntermd_2(i,j,k,iblock)
1083               do m=1,ntermd_2(i,j,k,iblock)
1084                 v2s(l,m,i,j,k,iblock)=v2s_all(l,m,i,j,k,iblock,iparm)
1085               enddo
1086             enddo
1087           enddo
1088         enddo
1089       enddo
1090       enddo
1091 ! Restore parameters of the cumulants
1092       do i=-nloctyp,nloctyp
1093         do j=1,2
1094           b1(j,i)=b1_all(j,i,iparm)
1095           b1tilde(j,i)=b1tilde_all(j,i,iparm)
1096           b2(j,i)=b2_all(j,i,iparm)
1097         enddo
1098         do j=1,2
1099           do k=1,2
1100             cc(k,j,i)=cc_all(k,j,i,iparm)
1101             ctilde(k,j,i)=ctilde_all(k,j,i,iparm)
1102             dd(k,j,i)=dd_all(k,j,i,iparm)
1103             dtilde(k,j,i)=dtilde_all(k,j,i,iparm)
1104             ee(k,j,i)=ee_all(k,j,i,iparm)
1105           enddo
1106         enddo
1107       enddo
1108 ! Restore the parameters of electrostatic interactions
1109       do i=1,2
1110         do j=1,2
1111           app(j,i)=app_all(j,i,iparm)
1112           bpp(j,i)=bpp_all(j,i,iparm)
1113           ael6(j,i)=ael6_all(j,i,iparm)
1114           ael3(j,i)=ael3_all(j,i,iparm)
1115         enddo
1116       enddo
1117 ! Restore sidechain parameters
1118       do i=1,ntyp
1119         do j=1,ntyp
1120           aa_aq(j,i)=aa_all(j,i,iparm)
1121           bb_aq(j,i)=bb_all(j,i,iparm)
1122           r0(j,i)=r0_all(j,i,iparm)
1123           sigma(j,i)=sigma_all(j,i,iparm)
1124           chi(j,i)=chi_all(j,i,iparm)
1125           augm(j,i)=augm_all(j,i,iparm)
1126           eps(j,i)=eps_all(j,i,iparm)
1127         enddo
1128       enddo
1129       do i=1,ntyp
1130         chip(i)=chip_all(i,iparm)
1131         alp(i)=alp_all(i,iparm)
1132       enddo
1133 ! Restore the SCp parameters
1134       do i=1,ntyp
1135         do j=1,2
1136           aad(i,j)=aad_all(i,j,iparm)
1137           bad(i,j)=bad_all(i,j,iparm)
1138         enddo
1139       enddo
1140 ! Restore disulfide-bond parameters
1141       ebr=ebr_all(iparm)
1142       d0cm=d0cm_all(iparm)
1143       akcm=akcm_all(iparm)
1144       akth=akth_all(iparm)
1145       akct=akct_all(iparm)
1146       v1ss=v1ss_all(iparm)
1147       v2ss=v2ss_all(iparm)
1148       v3ss=v3ss_all(iparm)
1149 ! Restore SC-backbone correlation parameters
1150       do i=-nsccortyp,nsccortyp
1151        do j=-nsccortyp,nsccortyp
1152
1153       nterm_sccor(j,i)=nterm_sccor_all(j,i,iparm)
1154         do l=1,3
1155            do k=1,nterm_sccor(j,i)
1156             v1sccor(k,l,j,i)=v1sccor_all(k,l,j,i,iparm)
1157             v2sccor(k,l,j,i)=v2sccor_all(k,l,j,i,iparm)
1158            enddo
1159           enddo
1160         enddo
1161       enddo
1162       return
1163       end subroutine restore_parm
1164 !--------------------------------------------------------------------------
1165 ! make_ensemble1.F
1166 !--------------------------------------------------------------------------
1167       subroutine make_ensembles(islice,*)
1168 ! construct the conformational ensembles at REMD temperatures
1169       use geometry_data, only:c
1170       use io_base, only:ilen
1171       use io_wham, only:pdboutW
1172 !      implicit none
1173 !      include "DIMENSIONS"
1174 !      include "DIMENSIONS.ZSCOPT"
1175 !      include "DIMENSIONS.FREE"
1176 #ifdef MPI
1177       include "mpif.h"
1178 !      include "COMMON.MPI"
1179       integer :: ierror,errcode,status(MPI_STATUS_SIZE)
1180 #endif
1181 !      include "COMMON.IOUNITS"
1182 !      include "COMMON.CONTROL"
1183 !      include "COMMON.FREE"
1184 !      include "COMMON.ENERGIES"
1185 !      include "COMMON.FFIELD"
1186 !      include "COMMON.INTERACT"
1187 !      include "COMMON.SBRIDGE"
1188 !      include "COMMON.CHAIN"
1189 !      include "COMMON.PROTFILES"
1190 !      include "COMMON.PROT"
1191       real(kind=4) :: csingle(3,nres*2)
1192       real(kind=8),dimension(6) :: fT,fTprim,fTbis
1193       real(kind=8) :: quot,quotl1,quotl,kfacl,&
1194         eprim,ebis,temper,kfac=2.4d0,T0=300.0d0
1195       real(kind=8) :: etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,&
1196             escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,&
1197             eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt, &
1198             ecation_prot, ecationcation
1199       integer :: i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
1200       real(kind=8) :: qfree,sumprob,eini,efree,rmsdev
1201       character(len=80) :: bxname
1202       character(len=2) :: licz1,licz2
1203       character(len=3) :: licz3,licz4
1204       character(len=5) :: ctemper
1205 !el      integer ilen
1206 !el      external ilen
1207       real(kind=4) :: Fdimless(MaxStr),Fdimless_(MaxStr)
1208       real(kind=8) :: enepot(MaxStr)
1209       integer :: iperm(MaxStr)
1210       integer :: islice
1211       integer,dimension(0:nprocs) :: scount_
1212
1213 #ifdef MPI
1214       if (me.eq.Master) then
1215 #endif
1216       write (licz2,'(bz,i2.2)') islice
1217       if (nslice.eq.1) then
1218         if (.not.separate_parset) then
1219           bxname = prefix(:ilen(prefix))//".bx"
1220         else
1221           write (licz3,'(bz,i3.3)') myparm
1222           bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx"
1223         endif
1224       else
1225         if (.not.separate_parset) then
1226           bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
1227         else
1228           write (licz3,'(bz,i3.3)') myparm
1229           bxname = prefix(:ilen(prefix))//"par_"//licz3// &
1230             "_slice_"//licz2//".bx"
1231         endif
1232       endif
1233       open (ientout,file=bxname,status="unknown",&
1234         form="unformatted",access="direct",recl=lenrec1)
1235 #ifdef MPI
1236       endif
1237 #endif
1238       do iparm=1,iparm
1239         if (iparm.ne.iparmprint) exit
1240         call restore_parm(iparm)
1241         do ib=1,nT_h(iparm)
1242 #ifdef DEBUG
1243           write (iout,*) "iparm",iparm," ib",ib
1244 #endif
1245           temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1246 !          quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1247 !          quotl=1.0d0
1248 !          kfacl=1.0d0
1249 !          do l=1,5
1250 !            quotl1=quotl
1251 !            quotl=quotl*quot
1252 !            kfacl=kfacl*kfac
1253 !            fT(l)=kfacl/(kfacl-1.0d0+quotl)
1254 !          enddo
1255 !el old rescale weights
1256 !
1257 !            if (rescale_mode.eq.1) then
1258 !              quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1259 #if defined(FUNCTH)
1260               tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1261               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
1262 #elif defined(FUNCT)
1263               ft(6)=quot
1264 #else
1265               ft(6)=1.0d0
1266 #endif
1267 !              quotl=1.0d0
1268 !              kfacl=1.0d0
1269 !              do l=1,5
1270 !                quotl1=quotl
1271 !                quotl=quotl*quot
1272 !                kfacl=kfacl*kfac
1273 !                fT(l)=kfacl/(kfacl-1.0d0+quotl)
1274 !              enddo
1275 !            else if (rescale_mode.eq.2) then
1276 !              quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1277 !#if defined(FUNCTH)
1278 !              tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1279 !              ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0
1280 !#elif defined(FUNCT)
1281 !              ft(6)=quot
1282 !#else
1283 !              ft(6)=1.0d0
1284 !#endif
1285 !              quotl=1.0d0
1286 !              do l=1,5
1287 !                quotl=quotl*quot
1288 !                fT(l)=1.12692801104297249644d0/ &
1289 !                   dlog(dexp(quotl)+dexp(-quotl))
1290 !              enddo
1291 !              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
1292 !            else if (rescale_mode.eq.0) then
1293 !              do l=1,5
1294 !                fT(l)=0.0d0
1295 !              enddo
1296 !            else
1297 !              write (iout,*) &
1298 !              "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
1299 !              call flush(iout)
1300 !              return 1
1301 !            endif
1302 ! el end old rescale weihgts
1303           call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3))
1304
1305 #ifdef MPI
1306           do i=1,scount(me1)
1307 #else
1308           do i=1,ntot(islice)
1309 #endif
1310             evdw=enetb(1,i,iparm)
1311 !            evdw_t=enetb(21,i,iparm)
1312             evdw_t=enetb(20,i,iparm)
1313 #ifdef SCP14
1314 !            evdw2_14=enetb(17,i,iparm)
1315             evdw2_14=enetb(18,i,iparm)
1316             evdw2=enetb(2,i,iparm)+evdw2_14
1317 #else
1318             evdw2=enetb(2,i,iparm)
1319             evdw2_14=0.0d0
1320 #endif
1321 #ifdef SPLITELE
1322             ees=enetb(3,i,iparm)
1323             evdw1=enetb(16,i,iparm)
1324 #else
1325             ees=enetb(3,i,iparm)
1326             evdw1=0.0d0
1327 #endif
1328             ecorr=enetb(4,i,iparm)
1329             ecorr5=enetb(5,i,iparm)
1330             ecorr6=enetb(6,i,iparm)
1331             eel_loc=enetb(7,i,iparm)
1332             eello_turn3=enetb(8,i,iparm)
1333             eello_turn4=enetb(9,i,iparm)
1334             eello_turn6=enetb(10,i,iparm)
1335             ebe=enetb(11,i,iparm)
1336             escloc=enetb(12,i,iparm)
1337             etors=enetb(13,i,iparm)
1338             etors_d=enetb(14,i,iparm)
1339             ehpb=enetb(15,i,iparm)
1340
1341             estr=enetb(17,i,iparm)
1342 !            estr=enetb(18,i,iparm)
1343 !            esccor=enetb(19,i,iparm)
1344             esccor=enetb(21,i,iparm)
1345 !            edihcnstr=enetb(20,i,iparm)
1346             edihcnstr=enetb(19,i,iparm)
1347 !#ifdef SPLITELE
1348 !            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1349 !            +wvdwpp*evdw1 &
1350 !            +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1351 !            +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1352 !            +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1353 !            +ft(2)*wturn3*eello_turn3 &
1354 !            +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc &
1355 !            +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1356 !            +wbond*estr
1357 !#else
1358 !            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1359 !            +ft(1)*welec*(ees+evdw1) &
1360 !            +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1361 !            +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1362 !            +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1363 !            +ft(2)*wturn3*eello_turn3 &
1364 !            +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1365 !            +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1366 !            +wbond*estr
1367 !#endif
1368
1369 #ifdef SPLITELE
1370             etot=wsc*evdw+wscp*evdw2+welec*ees &
1371             +wvdwpp*evdw1 &
1372             +wang*ebe+wtor*etors+wscloc*escloc &
1373             +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
1374             +wcorr6*ecorr6+wturn4*eello_turn4 &
1375             +wturn3*eello_turn3 &
1376             +wturn6*eello_turn6+wel_loc*eel_loc &
1377             +edihcnstr+wtor_d*etors_d+wsccor*esccor &
1378             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
1379 #else
1380             etot=wsc*evdw+wscp*evdw2 &
1381             +welec*(ees+evdw1) &
1382             +wang*ebe+wtor*etors+wscloc*escloc &
1383             +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
1384             +wcorr6*ecorr6+wturn4*eello_turn4 &
1385             +wturn3*eello_turn3 &
1386             +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
1387             +wtor_d*etors_d+wsccor*esccor &
1388             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
1389 #endif
1390
1391 #ifdef MPI
1392             Fdimless(i)= &
1393               beta_h(ib,iparm)*etot-entfac(i)
1394             potE(i,iparm)=etot
1395 #ifdef DEBUG
1396             write (iout,*) i,indstart(me)+i-1,ib,&
1397              1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),&
1398              -entfac(i),Fdimless(i)
1399 #endif
1400 #else
1401             Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
1402             potE(i,iparm)=etot
1403 #endif
1404           enddo   ! i
1405 #ifdef MPI
1406           do i=1,scount(me1)
1407             Fdimless_(i)=Fdimless(i)
1408           enddo
1409           call MPI_Gatherv(Fdimless_(1),scount(me),&
1410            MPI_REAL,Fdimless(1),&
1411            scount_(0),idispl(0),MPI_REAL,Master,&
1412            WHAM_COMM, IERROR)
1413 #ifdef DEBUG
1414           call MPI_Gatherv(potE(1,iparm),scount_(me),&
1415            MPI_DOUBLE_PRECISION,potE(1,iparm),&
1416            scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1417            WHAM_COMM, IERROR)
1418           call MPI_Gatherv(entfac(1),scount(me),&
1419            MPI_DOUBLE_PRECISION,entfac(1),&
1420            scount_(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1421            WHAM_COMM, IERROR)
1422 #endif
1423           if (me.eq.Master) then
1424 #ifdef DEBUG
1425           write (iout,*) "The FDIMLESS array before sorting"
1426           do i=1,ntot(islice)
1427             write (iout,*) i,fdimless(i)
1428           enddo
1429 #endif
1430 #endif
1431           do i=1,ntot(islice)
1432             iperm(i)=i
1433           enddo
1434           call mysort1(ntot(islice),Fdimless,iperm)
1435 #ifdef DEBUG
1436           write (iout,*) "The FDIMLESS array after sorting"
1437           do i=1,ntot(islice)
1438             write (iout,*) i,iperm(i),fdimless(i)
1439           enddo
1440 #endif
1441           qfree=0.0d0
1442           do i=1,ntot(islice)
1443             qfree=qfree+exp(-fdimless(i)+fdimless(1))
1444           enddo
1445 !          write (iout,*) "qfree",qfree
1446           nlist=1
1447           sumprob=0.0
1448           do i=1,min0(ntot(islice),ensembles)
1449             sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
1450 #ifdef DEBUG
1451             write (iout,*) i,ib,beta_h(ib,iparm),&
1452              1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),&
1453              potE(iperm(i),iparm),&
1454              -entfac(iperm(i)),fdimless(i),sumprob
1455 #endif
1456             if (sumprob.gt.0.99d0) goto 122
1457             nlist=nlist+1
1458           enddo
1459   122     continue
1460 #ifdef MPI
1461           endif
1462           call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,&
1463              IERROR)
1464           call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,&
1465              IERROR)
1466           do i=1,nlist
1467             ii=iperm(i)
1468             iproc=0
1469             do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
1470               iproc=iproc+1
1471             enddo
1472             if (iproc.ge.nprocs) then
1473               write (iout,*) "Fatal error: processor out of range",iproc
1474               call flush(iout)
1475               if (bxfile) then
1476                 close (ientout)
1477               else
1478                 close (ientout,status="delete")
1479               endif
1480               return 1
1481             endif
1482             ik=ii-indstart(iproc)+1
1483             if (iproc.ne.Master) then
1484               if (me.eq.iproc) then
1485 #ifdef DEBUG
1486                 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,&
1487                  " energy",potE(ik,iparm)
1488 #endif
1489                 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,&
1490                   Master,i,WHAM_COMM,IERROR)
1491               else if (me.eq.Master) then
1492                 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,&
1493                   WHAM_COMM,STATUS,IERROR)
1494               endif
1495             else if (me.eq.Master) then
1496               enepot(i)=potE(ik,iparm)
1497             endif
1498           enddo
1499 #else
1500           do i=1,nlist
1501             enepot(i)=potE(iperm(i),iparm)
1502           enddo
1503 #endif
1504 #ifdef MPI
1505           if (me.eq.Master) then
1506 #endif
1507           write(licz3,'(bz,i3.3)') iparm
1508           write(licz2,'(bz,i2.2)') islice
1509           if (temper.lt.100.0d0) then
1510             write(ctemper,'(f3.0)') temper
1511           else if (temper.lt.1000.0) then
1512             write (ctemper,'(f4.0)') temper
1513           else
1514             write (ctemper,'(f5.0)') temper
1515           endif
1516           if (nparmset.eq.1) then
1517             if (separate_parset) then
1518               write(licz4,'(bz,i3.3)') myparm
1519               pdbname=prefix(:ilen(prefix))//"_par"//licz4
1520             else
1521               pdbname=prefix(:ilen(prefix))
1522             endif
1523           else
1524             pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
1525           endif
1526           if (nslice.eq.1) then
1527             pdbname=pdbname(:ilen(pdbname))//"_T_"// &
1528               ctemper(:ilen(ctemper))//"pdb"
1529           else
1530             pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"// &
1531               ctemper(:ilen(ctemper))//"pdb"
1532           endif
1533           open(ipdb,file=pdbname)
1534           do i=1,nlist
1535             read (ientout,rec=iperm(i)) &
1536               ((csingle(l,k),l=1,3),k=1,nres),&
1537               ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1538               nss,(ihpb(k),jhpb(k),k=1,nss),&
1539               eini,efree,rmsdev,iscor
1540             do j=1,2*nres
1541               do k=1,3
1542                 c(k,j)=csingle(k,j)
1543               enddo
1544             enddo
1545             eini=fdimless(i)
1546             call pdboutW(iperm(i),temper,eini,enepot(i),efree,rmsdev)
1547           enddo
1548 #ifdef MPI
1549           endif
1550 #endif
1551         enddo     ! ib
1552       enddo       ! iparm
1553       if (bxfile) then
1554         close(ientout)
1555       else
1556         close(ientout,status="delete")
1557       endif
1558       do i=0,nprocs
1559         scount(i)=scount_(i)
1560       enddo
1561       return
1562       end subroutine make_ensembles
1563 !--------------------------------------------------------------------------
1564       subroutine mysort1(n, x, ipermut)
1565 !      implicit none
1566       integer :: i,j,imax,ipm,n
1567       real(kind=4) :: x(n)
1568       integer :: ipermut(n)
1569       real(kind=4) :: xtemp
1570       do i=1,n
1571         xtemp=x(i)
1572         imax=i
1573         do j=i+1,n
1574           if (x(j).lt.xtemp) then
1575             imax=j
1576             xtemp=x(j)
1577           endif
1578         enddo
1579         x(imax)=x(i)
1580         x(i)=xtemp
1581         ipm=ipermut(imax)
1582         ipermut(imax)=ipermut(i)
1583         ipermut(i)=ipm
1584       enddo
1585       return
1586       end subroutine mysort1
1587 !--------------------------------------------------------------------------
1588       subroutine alloc_enecalc_arrays(iparm)
1589
1590       use control_data
1591       use geometry_data, only:maxlob
1592       integer :: iparm
1593 !---------------------------
1594 ! COMMON.ENERGIES form wham_data
1595 !      common /energies/
1596       allocate(potE(MaxStr_Proc,iparm)) !(MaxStr_Proc,Max_Parm)
1597       allocate(entfac(MaxStr_Proc)) !(MaxStr_Proc)
1598       allocate(q(nQ+2,MaxStr_Proc)) !(MaxQ+2,MaxStr_Proc)
1599       allocate(enetb(max_ene,MaxStr_Proc,iparm)) !(max_ene,MaxStr_Proc,Max_Parm)
1600 !
1601 ! allocate ENECALC arrays
1602 !---------------------------
1603 ! COMMON.ALLPARM
1604 !      common /allparm/
1605       allocate(ww_all(max_eneW,iparm)) !(max_ene,max_parm) ! max_eneW
1606       allocate(vbldp0_all(iparm),akp_all(nParmSet)) !(max_parm)
1607       allocate(vbldsc0_all(maxbondterm,ntyp,iparm),&
1608         aksc_all(maxbondterm,ntyp,iparm),&
1609         abond0_all(maxbondterm,ntyp,iparm)) !(maxbondterm,ntyp,max_parm)
1610       allocate(a0thet_all(-ntyp:ntyp,iparm)) !(-ntyp:ntyp,max_parm)
1611       allocate(athet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm),&
1612         bthet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm)) !(2,-ntyp:ntyp,-1:1,-1:1,max_parm)
1613       allocate(polthet_all(0:3,-ntyp:ntyp,iparm)) !(0:3,-ntyp:ntyp,max_parm)
1614       allocate(gthet_all(3,-ntyp:ntyp,iparm)) !(3,-ntyp:ntyp,max_parm)
1615       allocate(theta0_all(-ntyp:ntyp,iparm),&
1616         sig0_all(-ntyp:ntyp,iparm),sigc0_all(-ntyp:ntyp,nParmSet)) !(-ntyp:ntyp,max_parm)
1617       allocate(aa0thet_all(-nthetyp-1:nthetyp+1,&
1618         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1619 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1620       allocate(aathet_all(ntheterm,-nthetyp-1:nthetyp+1,&
1621         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1622 !(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1623       allocate(bbthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1624         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1625       allocate(ccthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1626         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1627       allocate(ddthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1628         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1629       allocate(eethet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1630         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1631 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,
1632 !     & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1633       allocate(ffthet_all1(ndouble,ndouble,ntheterm3,&
1634         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1635         -nthetyp-1:nthetyp+1,iparm))
1636       allocate(ggthet_all1(ndouble,ndouble,ntheterm3,&
1637         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1638         -nthetyp-1:nthetyp+1,iparm))
1639       allocate(ffthet_all2(ndouble,ndouble,ntheterm3,&
1640         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1641         -nthetyp-1:nthetyp+1,iparm))
1642       allocate(ggthet_all2(ndouble,ndouble,ntheterm3,&
1643         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1644         -nthetyp-1:nthetyp+1,iparm))
1645 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1646 !-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm)
1647       allocate(dsc_all(ntyp1,iparm),dsc0_all(ntyp1,nParmSet)) !(ntyp1,max_parm)
1648       allocate(bsc_all(maxlob,ntyp,iparm))
1649 !(maxlob,ntyp,max_parm)
1650       allocate(censc_all(3,maxlob,-ntyp:ntyp,iparm)) !(3,maxlob,-ntyp:ntyp,max_parm)
1651       allocate(gaussc_all(3,3,maxlob,-ntyp:ntyp,iparm)) !(3,3,maxlob,-ntyp:ntyp,max_parm)
1652       allocate(sc_parmin_all(65,ntyp,iparm)) !(65,ntyp,max_parm)
1653       allocate(v0_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1654 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1655       allocate(v1_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1656       allocate(v2_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1657 !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1658       allocate(vlor1_all(maxlor,ntortyp,ntortyp,iparm))
1659       allocate(vlor2_all(maxlor,ntortyp,ntortyp,iparm))
1660       allocate(vlor3_all(maxlor,ntortyp,ntortyp,iparm)) !(maxlor,maxtor,maxtor,max_parm)
1661       allocate(v1c_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1662         -ntortyp:ntortyp,2,iparm))
1663       allocate(v1s_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1664         -ntortyp:ntortyp,2,iparm))
1665 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1666       allocate(v2c_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1667         -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1668 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1669       allocate(v2s_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1670         -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1671 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1672       allocate(b1_all(2,-nloctyp:nloctyp,iparm))
1673       allocate(b2_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1674       allocate(cc_all(2,2,-nloctyp:nloctyp,iparm))
1675       allocate(dd_all(2,2,-nloctyp:nloctyp,iparm))
1676       allocate(ee_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1677       allocate(ctilde_all(2,2,-nloctyp:nloctyp,iparm))
1678       allocate(dtilde_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1679       allocate(b1tilde_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1680       allocate(app_all(2,2,iparm),bpp_all(2,2,nParmSet),&
1681         ael6_all(2,2,iparm),ael3_all(2,2,nParmSet)) !(2,2,max_parm)
1682       allocate(aad_all(ntyp,2,iparm),bad_all(ntyp,2,nParmSet)) !(ntyp,2,max_parm)
1683       allocate(aa_all(ntyp,ntyp,iparm),bb_all(ntyp,ntyp,nParmSet),&
1684         augm_all(ntyp,ntyp,iparm),eps_all(ntyp,ntyp,nParmSet),&
1685         sigma_all(ntyp,ntyp,iparm),r0_all(ntyp,ntyp,nParmSet),&
1686         chi_all(ntyp,ntyp,iparm)) !(ntyp,ntyp,max_parm)
1687       allocate(chip_all(ntyp,iparm),alp_all(ntyp,nParmSet)) !(ntyp,max_parm)
1688       allocate(ebr_all(iparm),d0cm_all(nParmSet),akcm_all(nParmSet),&
1689         akth_all(iparm),akct_all(nParmSet),v1ss_all(nParmSet),&
1690         v2ss_all(iparm),v3ss_all(nParmSet)) !(max_parm)
1691       allocate(v1sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1692       allocate(v2sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1693 !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
1694       allocate(nlob_all(ntyp1,iparm)) !(ntyp1,max_parm)
1695       allocate(nlor_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1696       allocate(nterm_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1697 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1698       allocate(ntermd1_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1699         -ntortyp:ntortyp,2,iparm))
1700       allocate(ntermd2_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1701         -ntortyp:ntortyp,2,iparm))
1702 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1703       allocate(nbondterm_all(ntyp,iparm)) !(ntyp,max_parm)
1704       allocate(ithetyp_all(-ntyp1:ntyp1,iparm)) !(-ntyp1:ntyp1,max_parm)
1705       allocate(nthetyp_all(iparm),ntheterm_all(nParmSet),&
1706         ntheterm2_all(iparm),ntheterm3_all(nParmSet),&
1707         nsingle_all(iparm),&
1708         ndouble_all(iparm),nntheterm_all(nParmSet)) !(max_parm)
1709       allocate(nterm_sccor_all(-ntyp:ntyp,-ntyp:ntyp,iparm)) !(-ntyp:ntyp,-ntyp:ntyp,max_parm)
1710 !
1711       end subroutine alloc_enecalc_arrays
1712 !--------------------------------------------------------------------------
1713 !--------------------------------------------------------------------------
1714       end module ene_calc