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=100000
 137 !el      real(kind=8) :: rmsnat,gyrate
 138 !el      external rmsnat,gyrate
 139       real(kind=8) :: tole=0.0d0
 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,max_ene
 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             elipbond,elipang,eliplj,elipelec,ecat_prottran,ecation_protang, &
1236             eliptran,ehomology_constr
1237
1238       integer :: i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
1239       real(kind=8) :: qfree,sumprob,eini,efree,rmsdev
1240       character(len=80) :: bxname
1241       character(len=2) :: licz1,licz2
1242       character(len=3) :: licz3,licz4
1243       character(len=5) :: ctemper
1244 !el      integer ilen
1245 !el      external ilen
1246       real(kind=4) :: Fdimless(MaxStr),Fdimless_(MaxStr)
1247       real(kind=8) :: enepot(MaxStr)
1248       integer :: iperm(MaxStr)
1249       integer :: islice
1250       integer,dimension(0:nprocs) :: scount_
1251       write(iout,*) "Begin make ensemble"
1252 #ifdef MPI
1253       if (me.eq.Master) then
1254 #endif
1255       write (licz2,'(bz,i2.2)') islice
1256       if (nslice.eq.1) then
1257         if (.not.separate_parset) then
1258           bxname = prefix(:ilen(prefix))//".bx"
1259         else
1260           write (licz3,'(bz,i3.3)') myparm
1261           bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx"
1262         endif
1263       else
1264         if (.not.separate_parset) then
1265           bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
1266         else
1267           write (licz3,'(bz,i3.3)') myparm
1268           bxname = prefix(:ilen(prefix))//"par_"//licz3// &
1269             "_slice_"//licz2//".bx"
1270         endif
1271       endif
1272       open (ientout,file=bxname,status="unknown",&
1273         form="unformatted",access="direct",recl=lenrec1)
1274 #ifdef MPI
1275       endif
1276 #endif
1277       write(iout,*) "iparmprint",iparmprint,iparm
1278       do iparm=1,iparm
1279         if (iparm.ne.iparmprint) exit
1280         call restore_parm(iparm)
1281         do ib=1,nT_h(iparm)
1282 !#ifdef DEBUG
1283           write (iout,*) "iparm",iparm," ib",ib
1284 !#endif
1285           temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1286 !          quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1287 !          quotl=1.0d0
1288 !          kfacl=1.0d0
1289 !          do l=1,5
1290 !            quotl1=quotl
1291 !            quotl=quotl*quot
1292 !            kfacl=kfacl*kfac
1293 !            fT(l)=kfacl/(kfacl-1.0d0+quotl)
1294 !          enddo
1295 !el old rescale weights
1296 !
1297 !            if (rescale_mode.eq.1) then
1298 !              quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1299 #if defined(FUNCTH)
1300               tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1301               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
1302 #elif defined(FUNCT)
1303               ft(6)=quot
1304 #else
1305               ft(6)=1.0d0
1306 #endif
1307 !              quotl=1.0d0
1308 !              kfacl=1.0d0
1309 !              do l=1,5
1310 !                quotl1=quotl
1311 !                quotl=quotl*quot
1312 !                kfacl=kfacl*kfac
1313 !                fT(l)=kfacl/(kfacl-1.0d0+quotl)
1314 !              enddo
1315 !            else if (rescale_mode.eq.2) then
1316 !              quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1317 !#if defined(FUNCTH)
1318 !              tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1319 !              ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0
1320 !#elif defined(FUNCT)
1321 !              ft(6)=quot
1322 !#else
1323 !              ft(6)=1.0d0
1324 !#endif
1325 !              quotl=1.0d0
1326 !              do l=1,5
1327 !                quotl=quotl*quot
1328 !                fT(l)=1.12692801104297249644d0/ &
1329 !                   dlog(dexp(quotl)+dexp(-quotl))
1330 !              enddo
1331 !              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
1332 !            else if (rescale_mode.eq.0) then
1333 !              do l=1,5
1334 !                fT(l)=0.0d0
1335 !              enddo
1336 !            else
1337 !              write (iout,*) &
1338 !              "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
1339 !              call flush(iout)
1340 !              return 1
1341 !            endif
1342 ! el end old rescale weihgts
1343           call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3))
1344
1345 #ifdef MPI
1346           do i=1,scount(me1)
1347 #else
1348           do i=1,ntot(islice)
1349 #endif
1350             evdw=enetb(1,i,iparm)
1351 !            evdw_t=enetb(21,i,iparm)
1352             evdw_t=enetb(20,i,iparm)
1353 #ifdef SCP14
1354 !            evdw2_14=enetb(17,i,iparm)
1355             evdw2_14=enetb(18,i,iparm)
1356             evdw2=enetb(2,i,iparm)+evdw2_14
1357 #else
1358             evdw2=enetb(2,i,iparm)
1359             evdw2_14=0.0d0
1360 #endif
1361 #ifdef SPLITELE
1362             ees=enetb(3,i,iparm)
1363             evdw1=enetb(16,i,iparm)
1364 #else
1365             ees=enetb(3,i,iparm)
1366             evdw1=0.0d0
1367 #endif
1368             ecorr=enetb(4,i,iparm)
1369             ecorr5=enetb(5,i,iparm)
1370             ecorr6=enetb(6,i,iparm)
1371             eel_loc=enetb(7,i,iparm)
1372             eello_turn3=enetb(8,i,iparm)
1373             eello_turn4=enetb(9,i,iparm)
1374             eello_turn6=enetb(10,i,iparm)
1375             ebe=enetb(11,i,iparm)
1376             escloc=enetb(12,i,iparm)
1377             etors=enetb(13,i,iparm)
1378             etors_d=enetb(14,i,iparm)
1379             ehpb=enetb(15,i,iparm)
1380
1381             estr=enetb(17,i,iparm)
1382 !            estr=enetb(18,i,iparm)
1383 !            esccor=enetb(19,i,iparm)
1384             esccor=enetb(21,i,iparm)
1385 !            edihcnstr=enetb(20,i,iparm)
1386             edihcnstr=enetb(19,i,iparm)
1387             ecationcation=enetb(42,i,iparm)
1388             ecation_prot=enetb(41,i,iparm)
1389             evdwpp =      enetb(26,i,iparm)
1390             eespp =      enetb(27,i,iparm)
1391             evdwpsb =      enetb(28,i,iparm)
1392             eelpsb =      enetb(29,i,iparm)
1393             evdwsb =      enetb(30,i,iparm)
1394             eelsb =      enetb(31,i,iparm)
1395             estr_nucl =      enetb(32,i,iparm)
1396             ebe_nucl =      enetb(33,i,iparm)
1397             esbloc =      enetb(34,i,iparm)
1398             etors_nucl =      enetb(35,i,iparm)
1399             etors_d_nucl =      enetb(36,i,iparm)
1400             ecorr_nucl =      enetb(37,i,iparm)
1401             ecorr3_nucl =      enetb(38,i,iparm)
1402             epeppho=   enetb(49,i,iparm)
1403             escpho=    enetb(48,i,iparm)
1404             epepbase=  enetb(47,i,iparm)
1405             escbase=   enetb(46,i,iparm)
1406             ecation_nucl=enetb(50,i,iparm)
1407             elipbond=enetb(52,i,iparm)
1408             elipang=enetb(53,i,iparm)
1409             eliplj=enetb(54,i,iparm)
1410             elipelec=enetb(55,i,iparm)
1411             ecat_prottran=enetb(56,i,iparm)
1412             ecation_protang=enetb(57,i,iparm)
1413             eliptran=enetb(22,i,iparm)
1414             ehomology_constr=enetb(51,i,iparm)
1415           if (homol_nset.gt.1) &
1416             ehomology_constr=waga_homology(homol_nset)*ehomology_constr
1417 !      wscbase=ww(46)
1418 !      wpepbase=ww(47)
1419 !      wscpho=ww(48)
1420 !      wpeppho=ww(49)
1421
1422 !#ifdef SPLITELE
1423 !            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1424 !            +wvdwpp*evdw1 &
1425 !            +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1426 !            +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1427 !            +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1428 !            +ft(2)*wturn3*eello_turn3 &
1429 !            +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc &
1430 !            +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1431 !            +wbond*estr
1432 !#else
1433 !            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1434 !            +ft(1)*welec*(ees+evdw1) &
1435 !            +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1436 !            +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1437 !            +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1438 !            +ft(2)*wturn3*eello_turn3 &
1439 !            +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1440 !            +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1441 !            +wbond*estr
1442 !#endif
1443
1444 #ifdef SPLITELE
1445             etot=wsc*evdw+wscp*evdw2+welec*ees &
1446             +wvdwpp*evdw1 &
1447             +wang*ebe+wtor*etors+wscloc*escloc &
1448             +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
1449             +wcorr6*ecorr6+wturn4*eello_turn4 &
1450             +wturn3*eello_turn3 &
1451             +wturn6*eello_turn6+wel_loc*eel_loc &
1452             +edihcnstr+wtor_d*etors_d+wsccor*esccor &
1453             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
1454             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1455             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1456             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1457             +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
1458             +wscbase*escbase+ehomology_constr&
1459             +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
1460             +wcatnucl*ecation_nucl&
1461             +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang+&
1462              wliptran*eliptran
1463
1464
1465 #else
1466             etot=wsc*evdw+wscp*evdw2 &
1467             +welec*(ees+evdw1) &
1468             +wang*ebe+wtor*etors+wscloc*escloc &
1469             +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
1470             +wcorr6*ecorr6+wturn4*eello_turn4 &
1471             +wturn3*eello_turn3 &
1472             +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
1473             +wtor_d*etors_d+wsccor*esccor &
1474             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
1475             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1476             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1477             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1478             +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
1479             +wscbase*escbase+ehomology_constr&
1480             +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
1481             +wcatnucl*ecation_nucl&
1482             +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
1483             +wliptran*eliptran
1484
1485
1486
1487 #endif
1488
1489 #ifdef MPI
1490             Fdimless(i)= &
1491               beta_h(ib,iparm)*etot-entfac(i)
1492             potE(i,iparm)=etot
1493 !#ifdef DEBUG
1494             write (iout,*) i,indstart(me)+i-1,ib,&
1495              1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),&
1496              -entfac(i),Fdimless(i)
1497 !#endif
1498 #else
1499             Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
1500             potE(i,iparm)=etot
1501 #endif
1502           enddo   ! i
1503 #ifdef MPI
1504           scount_(:)=0
1505           do i=1,scount(me1)
1506             Fdimless_(i)=Fdimless(i)
1507           enddo
1508           write(iout,*) "before gather Fdimless"
1509           write(iout,*) scount(me),scount_(0),idispl(0)
1510           write (iout,*) "added update of scount_"
1511          call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount_(0), 1, &
1512         MPI_INTEGER, WHAM_COMM, IERROR)
1513
1514
1515           call MPI_Gatherv(Fdimless_(1),scount(me),&
1516            MPI_REAL,Fdimless(1),&
1517            scount_(0),idispl(0),MPI_REAL,Master,&
1518            WHAM_COMM, IERROR)
1519           write(iout,*) "after gather Fdimless"
1520 #ifdef DEBUG
1521           call MPI_Gatherv(potE(1,iparm),scount_(me),&
1522            MPI_DOUBLE_PRECISION,potE(1,iparm),&
1523            scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1524            WHAM_COMM, IERROR)
1525           call MPI_Gatherv(entfac(1),scount(me),&
1526            MPI_DOUBLE_PRECISION,entfac(1),&
1527            scount_(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1528            WHAM_COMM, IERROR)
1529 #endif
1530           if (me.eq.Master) then
1531 !#ifdef DEBUG
1532           write (iout,*) "The FDIMLESS array before sorting"
1533           do i=1,ntot(islice)
1534             write (iout,*) i,fdimless(i)
1535           enddo
1536 !#endif
1537 #endif
1538           do i=1,ntot(islice)
1539             iperm(i)=i
1540           enddo
1541           call mysort1(ntot(islice),Fdimless,iperm)
1542 !#ifdef DEBUG
1543           write (iout,*) "The FDIMLESS array after sorting"
1544           do i=1,ntot(islice)
1545             write (iout,*) i,iperm(i),fdimless(i)
1546           enddo
1547 !#endif
1548           qfree=0.0d0
1549           do i=1,ntot(islice)
1550             qfree=qfree+exp(-fdimless(i)+fdimless(1))
1551           enddo
1552 !          write (iout,*) "qfree",qfree
1553           nlist=1
1554           sumprob=0.0
1555           do i=1,min0(ntot(islice),ensembles)
1556             sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
1557 !#ifdef DEBUG
1558             write (iout,*) i,ib,beta_h(ib,iparm),&
1559              1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),&
1560              potE(iperm(i),iparm),&
1561              -entfac(iperm(i)),fdimless(i),sumprob
1562 !#endif
1563             if (sumprob.gt.0.99d0) goto 122
1564             nlist=nlist+1
1565           enddo
1566   122     continue
1567 #ifdef MPI
1568           endif
1569           call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,&
1570              IERROR)
1571           call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,&
1572              IERROR)
1573           do i=1,nlist
1574             ii=iperm(i)
1575             iproc=0
1576             do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
1577               iproc=iproc+1
1578             enddo
1579             if (iproc.ge.nprocs) then
1580               write (iout,*) "Fatal error: processor out of range",iproc
1581               call flush(iout)
1582               if (bxfile) then
1583                 close (ientout)
1584               else
1585                 close (ientout,status="delete")
1586               endif
1587               return 1
1588             endif
1589             ik=ii-indstart(iproc)+1
1590             if (iproc.ne.Master) then
1591               if (me.eq.iproc) then
1592 !#ifdef DEBUG
1593                 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,&
1594                  " energy",potE(ik,iparm)
1595 !#endif
1596                 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,&
1597                   Master,i,WHAM_COMM,IERROR)
1598               else if (me.eq.Master) then
1599                 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,&
1600                   WHAM_COMM,STATUS,IERROR)
1601               endif
1602             else if (me.eq.Master) then
1603               enepot(i)=potE(ik,iparm)
1604             endif
1605           enddo
1606 #else
1607           do i=1,nlist
1608             enepot(i)=potE(iperm(i),iparm)
1609           enddo
1610 #endif
1611           write(iout,*) "DEBUG",me
1612 #ifdef MPI
1613           if (me.eq.Master) then
1614 #endif
1615           write(licz3,'(bz,i3.3)') iparm
1616           write(licz2,'(bz,i2.2)') islice
1617           if (temper.lt.100.0d0) then
1618             write(ctemper,'(f3.0)') temper
1619           else if (temper.lt.1000.0) then
1620             write (ctemper,'(f4.0)') temper
1621           else
1622             write (ctemper,'(f5.0)') temper
1623           endif
1624           if (nparmset.eq.1) then
1625             if (separate_parset) then
1626               write(licz4,'(bz,i3.3)') myparm
1627               pdbname=prefix(:ilen(prefix))//"_par"//licz4
1628             else
1629               pdbname=prefix(:ilen(prefix))
1630             endif
1631           else
1632             pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
1633           endif
1634           if (nslice.eq.1) then
1635             pdbname=pdbname(:ilen(pdbname))//"_T_"// &
1636               ctemper(:ilen(ctemper))//"pdb"
1637           else
1638             pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"// &
1639               ctemper(:ilen(ctemper))//"pdb"
1640           endif
1641           open(ipdb,file=pdbname)
1642           do i=1,nlist
1643             read (ientout,rec=iperm(i)) &
1644               ((csingle(l,k),l=1,3),k=1,nres),&
1645               ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1646               nss,(ihpb(k),jhpb(k),k=1,nss),&
1647               eini,efree,rmsdev,iscor
1648             do j=1,2*nres
1649               do k=1,3
1650                 c(k,j)=csingle(k,j)
1651               enddo
1652             enddo
1653             call returnbox
1654             do j=1,2*nres
1655               do k=1,3
1656                 csingle(k,j)=c(k,j)
1657               enddo
1658             enddo
1659
1660             eini=fdimless(i)
1661             call pdboutW(iperm(i),temper,eini,enepot(i),efree,rmsdev)
1662           enddo
1663 #ifdef MPI
1664           endif
1665 #endif
1666         enddo     ! ib
1667       enddo       ! iparm
1668       if (bxfile) then
1669         close(ientout)
1670       else
1671         close(ientout,status="delete")
1672       endif
1673       do i=0,nprocs
1674         scount(i)=scount_(i)
1675       enddo
1676       return
1677       end subroutine make_ensembles
1678 !--------------------------------------------------------------------------
1679       subroutine mysort1(n, x, ipermut)
1680 !      implicit none
1681       integer :: i,j,imax,ipm,n
1682       real(kind=4) :: x(n)
1683       integer :: ipermut(n)
1684       real(kind=4) :: xtemp
1685       do i=1,n
1686         xtemp=x(i)
1687         imax=i
1688         do j=i+1,n
1689           if (x(j).lt.xtemp) then
1690             imax=j
1691             xtemp=x(j)
1692           endif
1693         enddo
1694         x(imax)=x(i)
1695         x(i)=xtemp
1696         ipm=ipermut(imax)
1697         ipermut(imax)=ipermut(i)
1698         ipermut(i)=ipm
1699       enddo
1700       return
1701       end subroutine mysort1
1702 !--------------------------------------------------------------------------
1703       subroutine alloc_enecalc_arrays(iparm)
1704
1705       use control_data
1706       use geometry_data, only:maxlob
1707       integer :: iparm
1708 !---------------------------
1709 ! COMMON.ENERGIES form wham_data
1710 !      common /energies/
1711       allocate(potE(MaxStr_Proc,iparm)) !(MaxStr_Proc,Max_Parm)
1712       allocate(entfac(MaxStr_Proc)) !(MaxStr_Proc)
1713       allocate(q(nQ+2,MaxStr_Proc)) !(MaxQ+2,MaxStr_Proc)
1714       allocate(enetb(max_ene,MaxStr_Proc,iparm)) !(max_ene,MaxStr_Proc,Max_Parm)
1715 !
1716 ! allocate ENECALC arrays
1717 !---------------------------
1718 ! COMMON.ALLPARM
1719 !      common /allparm/
1720       allocate(ww_all(max_eneW,iparm)) !(max_ene,max_parm) ! max_eneW
1721       allocate(vbldp0_all(iparm),akp_all(nParmSet)) !(max_parm)
1722       allocate(vbldsc0_all(maxbondterm,ntyp,iparm),&
1723         aksc_all(maxbondterm,ntyp,iparm),&
1724         abond0_all(maxbondterm,ntyp,iparm)) !(maxbondterm,ntyp,max_parm)
1725       allocate(a0thet_all(-ntyp:ntyp,iparm)) !(-ntyp:ntyp,max_parm)
1726       allocate(athet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm),&
1727         bthet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm)) !(2,-ntyp:ntyp,-1:1,-1:1,max_parm)
1728       allocate(polthet_all(0:3,-ntyp:ntyp,iparm)) !(0:3,-ntyp:ntyp,max_parm)
1729       allocate(gthet_all(3,-ntyp:ntyp,iparm)) !(3,-ntyp:ntyp,max_parm)
1730       allocate(theta0_all(-ntyp:ntyp,iparm),&
1731         sig0_all(-ntyp:ntyp,iparm),sigc0_all(-ntyp:ntyp,nParmSet)) !(-ntyp:ntyp,max_parm)
1732       allocate(aa0thet_all(-nthetyp-1:nthetyp+1,&
1733         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1734 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1735       allocate(aathet_all(ntheterm,-nthetyp-1:nthetyp+1,&
1736         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1737 !(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1738       allocate(bbthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1739         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1740       allocate(ccthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1741         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1742       allocate(ddthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1743         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1744       allocate(eethet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1745         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1746 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,
1747 !     & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1748       allocate(ffthet_all1(ndouble,ndouble,ntheterm3,&
1749         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1750         -nthetyp-1:nthetyp+1,iparm))
1751       allocate(ggthet_all1(ndouble,ndouble,ntheterm3,&
1752         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1753         -nthetyp-1:nthetyp+1,iparm))
1754       allocate(ffthet_all2(ndouble,ndouble,ntheterm3,&
1755         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1756         -nthetyp-1:nthetyp+1,iparm))
1757       allocate(ggthet_all2(ndouble,ndouble,ntheterm3,&
1758         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1759         -nthetyp-1:nthetyp+1,iparm))
1760 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1761 !-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm)
1762       allocate(dsc_all(ntyp1,iparm),dsc0_all(ntyp1,nParmSet)) !(ntyp1,max_parm)
1763       allocate(bsc_all(maxlob,ntyp,iparm))
1764 !(maxlob,ntyp,max_parm)
1765       allocate(censc_all(3,maxlob,-ntyp:ntyp,iparm)) !(3,maxlob,-ntyp:ntyp,max_parm)
1766       allocate(gaussc_all(3,3,maxlob,-ntyp:ntyp,iparm)) !(3,3,maxlob,-ntyp:ntyp,max_parm)
1767       allocate(sc_parmin_all(65,ntyp,iparm)) !(65,ntyp,max_parm)
1768       allocate(v0_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1769 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1770       allocate(v1_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1771       allocate(v2_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1772 !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1773       allocate(vlor1_all(maxlor,ntortyp,ntortyp,iparm))
1774       allocate(vlor2_all(maxlor,ntortyp,ntortyp,iparm))
1775       allocate(vlor3_all(maxlor,ntortyp,ntortyp,iparm)) !(maxlor,maxtor,maxtor,max_parm)
1776       allocate(v1c_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1777         -ntortyp:ntortyp,2,iparm))
1778       allocate(v1s_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1779         -ntortyp:ntortyp,2,iparm))
1780 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1781       allocate(v2c_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1782         -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1783 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1784       allocate(v2s_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1785         -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1786 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1787       allocate(b1_all(2,-nloctyp:nloctyp,iparm))
1788       allocate(b2_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1789       allocate(cc_all(2,2,-nloctyp:nloctyp,iparm))
1790       allocate(dd_all(2,2,-nloctyp:nloctyp,iparm))
1791       allocate(ee_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1792       allocate(ctilde_all(2,2,-nloctyp:nloctyp,iparm))
1793       allocate(dtilde_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1794       allocate(b1tilde_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1795       allocate(app_all(2,2,iparm),bpp_all(2,2,nParmSet),&
1796         ael6_all(2,2,iparm),ael3_all(2,2,nParmSet)) !(2,2,max_parm)
1797       allocate(aad_all(ntyp,2,iparm),bad_all(ntyp,2,nParmSet)) !(ntyp,2,max_parm)
1798       allocate(aa_all(ntyp,ntyp,iparm),bb_all(ntyp,ntyp,nParmSet),&
1799         augm_all(ntyp,ntyp,iparm),eps_all(ntyp,ntyp,nParmSet),&
1800         sigma_all(ntyp,ntyp,iparm),r0_all(ntyp,ntyp,nParmSet),&
1801         chi_all(ntyp,ntyp,iparm)) !(ntyp,ntyp,max_parm)
1802       allocate(chip_all(ntyp,iparm),alp_all(ntyp,nParmSet)) !(ntyp,max_parm)
1803       allocate(ebr_all(iparm),d0cm_all(nParmSet),akcm_all(nParmSet),&
1804         akth_all(iparm),akct_all(nParmSet),v1ss_all(nParmSet),&
1805         v2ss_all(iparm),v3ss_all(nParmSet)) !(max_parm)
1806       allocate(v1sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1807       allocate(v2sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1808 !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
1809       allocate(nlob_all(ntyp1,iparm)) !(ntyp1,max_parm)
1810       allocate(nlor_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1811       allocate(nterm_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1812 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1813       allocate(ntermd1_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1814         -ntortyp:ntortyp,2,iparm))
1815       allocate(ntermd2_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1816         -ntortyp:ntortyp,2,iparm))
1817 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1818       allocate(nbondterm_all(ntyp,iparm)) !(ntyp,max_parm)
1819       allocate(ithetyp_all(-ntyp1:ntyp1,iparm)) !(-ntyp1:ntyp1,max_parm)
1820       allocate(nthetyp_all(iparm),ntheterm_all(nParmSet),&
1821         ntheterm2_all(iparm),ntheterm3_all(nParmSet),&
1822         nsingle_all(iparm),&
1823         ndouble_all(iparm),nntheterm_all(nParmSet)) !(max_parm)
1824       allocate(nterm_sccor_all(-ntyp:ntyp,-ntyp:ntyp,iparm)) !(-ntyp:ntyp,-ntyp:ntyp,max_parm)
1825 !
1826       end subroutine alloc_enecalc_arrays
1827 !--------------------------------------------------------------------------
1828 !--------------------------------------------------------------------------
1829       end module ene_calc