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