source/wham/enecalc.F90

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