source/wham/enecalc.F90

   1       module ene_calc
   2 !-----------------------------------------------------------------------------
   3       use io_units
   4       use wham_data
   5       use control_data, only: tor_mode
   6 !
   7       use geometry_data, only:nres,boxxsize,boxysize,boxzsize
   8       use energy_data
   9 !      use control_data, only:maxthetyp1
  10       use energy, only:etotal,enerprint,rescale_weights
  11 #ifdef MPI
  12       use MPI_data
  13 !      include "mpif.h"
  14 !      include "COMMON.MPI"
  15 #endif
  16       implicit none
  17 !-----------------------------------------------------------------------------
  18 ! COMMON.ALLPARM
  19 !      common /allparm/
  20       real(kind=8),dimension(:,:),allocatable :: ww_all !(max_ene,max_parm) ! max_eneW
  21       real(kind=8),dimension(:),allocatable :: vbldp0_all,akp_all !(max_parm)
  22       real(kind=8),dimension(:,:,:),allocatable :: vbldsc0_all,&
  23         aksc_all,abond0_all !(maxbondterm,ntyp,max_parm)
  24       real(kind=8),dimension(:,:),allocatable :: a0thet_all !(-ntyp:ntyp,max_parm)
  25       real(kind=8),dimension(:,:,:,:,:),allocatable :: athet_all,&
  26         bthet_all !(2,-ntyp:ntyp,-1:1,-1:1,max_parm)
  27       real(kind=8),dimension(:,:,:),allocatable :: polthet_all !(0:3,-ntyp:ntyp,max_parm)
  28       real(kind=8),dimension(:,:,:),allocatable :: gthet_all !(3,-ntyp:ntyp,max_parm)
  29       real(kind=8),dimension(:,:),allocatable :: theta0_all,&
  30         sig0_all,sigc0_all !(-ntyp:ntyp,max_parm)
  31       real(kind=8),dimension(:,:,:,:,:),allocatable :: aa0thet_all
  32 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
  33       real(kind=8),dimension(:,:,:,:,:,:),allocatable :: aathet_all
  34 !(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
  35       real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: bbthet_all,&
  36         ccthet_all,ddthet_all,eethet_all !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,
  37 !     & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
  38       real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: ffthet_all1,&
  39         ggthet_all1,ffthet_all2,ggthet_all2 !(maxdouble,maxdouble,maxtheterm3,
  40 !     &  -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm)
  41       real(kind=8),dimension(:,:),allocatable :: dsc_all,dsc0_all !(ntyp1,max_parm)
  42       real(kind=8),dimension(:,:,:),allocatable :: bsc_all !(maxlob,ntyp,max_parm)
  43       real(kind=8),dimension(:,:,:,:),allocatable :: censc_all !(3,maxlob,-ntyp:ntyp,max_parm)
  44       real(kind=8),dimension(:,:,:,:,:),allocatable :: gaussc_all !(3,3,maxlob,-ntyp:ntyp,max_parm)
  45       real(kind=8),dimension(:,:,:),allocatable :: sc_parmin_all !(65,ntyp,max_parm)
  46       real(kind=8),dimension(:,:,:,:),allocatable :: v0_all
  47 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
  48       real(kind=8),dimension(:,:,:,:,:),allocatable :: v1_all,&
  49         v2_all !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
  50       real(kind=8),dimension(:,:,:,:),allocatable :: vlor1_all,&
  51         vlor2_all,vlor3_all !(maxlor,maxtor,maxtor,max_parm)
  52       real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v1c_all,&
  53         v1s_all !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
  54       real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v2c_all
  55 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
  56       real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v2s_all
  57 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
  58       real(kind=8),dimension(:,:,:),allocatable :: b1_all,b2_all !(2,-maxtor:maxtor,max_parm)
  59       real(kind=8),dimension(:,:,:,:),allocatable :: cc_all,dd_all,&
  60         ee_all !(2,2,-maxtor:maxtor,max_parm)
  61       real(kind=8),dimension(:,:,:,:),allocatable :: ctilde_all,&
  62         dtilde_all !(2,2,-maxtor:maxtor,max_parm)
  63       real(kind=8),dimension(:,:,:),allocatable :: b1tilde_all !(2,-maxtor:maxtor,max_parm)
  64       real(kind=8),dimension(:,:,:),allocatable :: app_all,bpp_all,&
  65         ael6_all,ael3_all !(2,2,max_parm)
  66       real(kind=8),dimension(:,:,:),allocatable :: aad_all,&
  67         bad_all !(ntyp,2,max_parm)
  68       real(kind=8),dimension(:,:,:),allocatable :: aa_all,bb_all,&
  69         augm_all,eps_all,sigma_all,r0_all,chi_all !(ntyp,ntyp,max_parm)
  70       real(kind=8),dimension(:,:),allocatable :: chip_all,alp_all !(ntyp,max_parm)
  71       real(kind=8),dimension(:),allocatable :: ebr_all,d0cm_all,&
  72         akcm_all,akth_all,akct_all,v1ss_all,v2ss_all,v3ss_all !(max_parm)
  73       real(kind=8),dimension(:,:,:,:,:),allocatable :: v1sccor_all,&
  74         v2sccor_all !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
  75       integer,dimension(:,:),allocatable :: nlob_all !(ntyp1,max_parm)
  76       integer,dimension(:,:,:,:),allocatable :: nlor_all,&
  77         nterm_all !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
  78       integer,dimension(:,:,:,:,:),allocatable :: ntermd1_all,&
  79         ntermd2_all !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
  80       integer,dimension(:,:),allocatable :: nbondterm_all !(ntyp,max_parm)
  81       integer,dimension(:,:),allocatable :: ithetyp_all !(-ntyp1:ntyp1,max_parm)
  82       integer,dimension(:),allocatable :: nthetyp_all,ntheterm_all,&
  83         ntheterm2_all,ntheterm3_all,nsingle_all,ndouble_all,&
  84         nntheterm_all !(max_parm)
  85       integer,dimension(:,:,:),allocatable :: nterm_sccor_all !(-ntyp:ntyp,-ntyp:ntyp,max_parm)
  86 !-----------------------------------------------------------------------------
  87 !
  88 !
  89 !-----------------------------------------------------------------------------
  90       contains
  91 !-----------------------------------------------------------------------------
  92       subroutine enecalc(islice,*)
  93
  94       use names
  95       use control_data, only:indpdb
  96       use geometry_data, only:c,phi,theta,alph,omeg,deg2rad,anatemp,&
  97                               vbld,rad2deg,dc_norm,dc,vbld_inv
  98       use io_base, only:gyrate!,briefout
  99       use geometry, only:int_from_cart1,returnbox
 100       use io_wham, only:pdboutW
 101       use io_database, only:opentmp
 102       use conform_compar, only:qwolynes,rmsnat
 103 !      include "DIMENSIONS"
 104 !      include "DIMENSIONS.ZSCOPT"
 105 !      include "DIMENSIONS.FREE"
 106 #ifdef MPI
 107 !      use MPI_data
 108       include "mpif.h"
 109 !      include "COMMON.MPI"
 110 #endif
 111 !      include "COMMON.CHAIN"
 112 !      include "COMMON.IOUNITS"
 113 !      include "COMMON.PROTFILES"
 114 !      include "COMMON.NAMES"
 115 !      include "COMMON.VAR"
 116 !      include "COMMON.SBRIDGE"
 117 !      include "COMMON.GEO"
 118 !      include "COMMON.FFIELD"
 119 !      include "COMMON.ENEPS"
 120 !      include "COMMON.LOCAL"
 121 !      include "COMMON.WEIGHTS"
 122 !      include "COMMON.INTERACT"
 123 !      include "COMMON.FREE"
 124 !      include "COMMON.ENERGIES"
 125 !      include "COMMON.CONTROL"
 126 !      include "COMMON.TORCNSTR"
 127 !      implicit none
 128 #ifdef MPI
 129       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 130 #endif
 131       character(len=64) :: nazwa
 132       character(len=80) :: bxname
 133       character(len=3) :: liczba
 134 !el      real(kind=8) :: qwolynes
 135 !el      external qwolynes
 136       integer :: errmsg_count,maxerrmsg_count=100
 137 !el      real(kind=8) :: rmsnat,gyrate
 138 !el      external rmsnat,gyrate
 139       real(kind=8) :: tole=1.0d-1
 140       integer i,itj,ii,iii,j,k,l,licz
 141       integer ir,ib,ipar,iparm
 142       integer iscor,islice
 143       real(kind=4) :: csingle(3,nres*2)
 144       real(kind=8) :: energ
 145       real(kind=8) :: temp
 146 !el      integer ilen,iroof
 147 !el      external ilen,iroof
 148       real(kind=8) :: energia(0:n_ene),rmsdev,efree,eini
 149 !el      real(kind=8) :: energia(0:max_ene),rmsdev,efree,eini
 150       real(kind=8) :: fT(6),quot,quotl,kfacl,kfac=2.4d0,T0=3.0d2
 151       real(kind=8) :: tt
 152       integer :: snk_p(MaxR,MaxT_h,nParmSet)!Max_parm)
 153       logical :: lerr
 154       character(len=64) :: bprotfile_temp
 155
 156 !      integer :: rec
 157       integer,dimension(0:nprocs) :: scount_
 158 !el      real(kind=8) :: rmsnat
 159
 160       rescale_mode=rescale_modeW
 161
 162       call opentmp(islice,ientout,bprotfile_temp)
 163       iii=0
 164       ii=0
 165 !el
 166 !      iparm=1
 167       errmsg_count=0
 168       write (iout,*) "enecalc: nparmset ",nparmset
 169 #ifdef MPI
 170       do iparm=1,nParmSet
 171         do ib=1,nT_h(iparm)
 172           do i=1,nR(ib,iparm)
 173             snk_p(i,ib,iparm)=0
 174           enddo
 175         enddo
 176       enddo
 177       do i=indstart(me1),indend(me1)
 178 write(iout,*)"enecalc_ i indstart",i,indstart(me1),indend(me1)
 179 #else
 180       do iparm=1,nParmSet
 181         do ib=1,nT_h(iparm)
 182           do i=1,nR(ib,iparm)
 183             snk(i,ib,iparm)=0
 184           enddo
 185         enddo
 186       enddo
 187       do i=1,ntot
 188 write(iout,*)"enecalc_ i ntot",i,ntot
 189 #endif
 190         read(ientout,rec=i,err=101) &
 191           ((csingle(l,k),l=1,3),k=1,nres),&
 192           ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
 193           nss,(ihpb(k),jhpb(k),k=1,nss),&
 194           eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
 195 !el debug
 196 !write(iout,*)"co wczytuje"
 197 !          write(iout,*)((csingle(l,k),l=1,3),k=1,nres),&
 198 !          ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
 199 !          nss,(ihpb(k),jhpb(k),k=1,nss),&
 200 !          eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
 201 !el --------
 202
 203 !write(iout,*)"ipar",ib,ipar,1.0d0/(beta_h(ib,ipar)*1.987D-3)
 204          if (indpdb.gt.0) then
 205            do k=1,nres
 206              do l=1,3
 207                c(l,k)=csingle(l,k)
 208              enddo
 209            enddo
 210            do k=nnt,nct
 211              do l=1,3
 212                c(l,k+nres)=csingle(l,k+nres)
 213              enddo
 214            enddo
 215            call returnbox
 216            do k=1,nres
 217              do l=1,3
 218                csingle(l,k)=c(l,k)
 219              enddo
 220            enddo
 221            do k=nnt,nct
 222              do l=1,3
 223                csingle(l,k+nres)=c(l,k+nres)
 224              enddo
 225            enddo
 226            anatemp= 1.0d0/(beta_h(ib,ipar)*1.987D-3)
 227            q(nQ+1,iii+1)=rmsnat(iii+1)
 228          endif
 229          q(nQ+2,iii+1)=gyrate(iii+1)
 230 ! write(iout,*)"wczyt",anatemp,q(nQ+2,iii+1) !el
 231 !        fT=T0*beta_h(ib,ipar)*1.987D-3
 232 !        ft=2.0d0/(1.0d0+1.0d0/(T0*beta_h(ib,ipar)*1.987D-3))
 233 ! EL start old rescale
 234 !        if (rescale_modeW.eq.1) then
 235 !          quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
 236 !#if defined(FUNCTH)
 237 !          tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
 238 !          ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
 239 !#elif defined(FUNCT)
 240 !          ft(6)=quot
 241 !#else
 242 !          ft(6)=1.0d0
 243 !#endif
 244 !          quotl=1.0d0
 245 !          kfacl=1.0d0
 246 !          do l=1,5
 247 !            quotl=quotl*quot
 248 !            kfacl=kfacl*kfac
 249 !            fT(l)=kfacl/(kfacl-1.0d0+quotl)
 250 !          enddo
 251 !        else if (rescale_modeW.eq.2) then
 252 !          quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
 253 !#if defined(FUNCTH)
 254 !          tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
 255 !          ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
 256 !#elif defined(FUNCT)
 257 !          ft(6)=quot
 258 !#else
 259 !          ft(6)=1.0d0
 260 !#endif
 261 !          quotl=1.0d0
 262 !          do l=1,5
 263 !            quotl=quotl*quot
 264 !            fT(l)=1.12692801104297249644d0/ &
 265 !               dlog(dexp(quotl)+dexp(-quotl))
 266 !          enddo
 267 !        else if (rescale_modeW.eq.0) then
 268 !          do l=1,5
 269 !            fT(l)=1.0d0
 270 !          enddo
 271 !        else
 272 !          write (iout,*) "Error in ECECALC: wrong RESCALE_MODE",&
 273 !           rescale_modeW
 274 !          call flush(iout)
 275 !          return 1
 276 !        endif
 277 !EL end old rescele
 278 !        write (iout,*) "T",1.0d0/(beta_h(ib,ipar)*1.987D-3)," T0",T0,
 279 !     &   " kfac",kfac,"quot",quot," fT",fT
 280 #ifdef DEBUG
 281             write(iout,*)"weights"
 282             write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
 283             wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
 284             wtor_d,wsccor,wbond
 285 #endif
 286
 287         do j=1,2*nres
 288           do k=1,3
 289             c(k,j)=csingle(k,j)
 290           enddo
 291         enddo
 292         call int_from_cart1(.false.)
 293         ii=ii+1
 294
 295 !        call rescale_weights(1.0d0/(beta_h(ib,ipar)*1.987D-3))
 296        do iparm=1,nparmset
 297 #ifdef DEBUG
 298             write (iout,*) "before restore w=",1.0d0/(beta_h(ib,ipar)*1.987D-3)
 299             write(iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
 300             wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
 301             wtor_d,wsccor,wbond
 302 #endif
 303         call restore_parm(iparm)
 304         call rescale_weights(1.0d0/(beta_h(ib,ipar)*1.987D-3))
 305 #ifdef DEBUG
 306             write (iout,*) "before etot w=",1.0d0/(beta_h(ib,ipar)*1.987D-3)
 307             write(iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
 308             wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
 309             wtor_d,wsccor,wbond
 310 #endif
 311 !        call etotal(energia(0),fT)
 312         call etotal(energia(0))
 313 !write(iout,*)"check c and dc after etotal",1.0d0/(0.001987*beta_h(ib,ipar))
 314 !do k=1,2*nres+2
 315 !write(iout,*)k,"c=",(c(l,k),l=1,3)
 316 !write(iout,*)k,"dc=",(dc(l,k),l=1,3)
 317 !write(iout,*)k,"dc_norm=",(dc_norm(l,k),l=1,3)
 318 !enddo
 319 !do k=1,nres*2
 320 !write(iout,*)k,"vbld=",vbld(k)
 321 !write(iout,*)k,"vbld_inv=",vbld_inv(k)
 322 !enddo
 323
 324 !write(iout,*)"energia",(energia(j),j=0,n_ene)
 325 !write(iout,*)"enerprint tuz po call etotal"
 326         call enerprint(energia(0))
 327 #ifdef DEBUG
 328         write (iout,*) "Conformation",i
 329           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 330           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 331 !        call enerprint(energia(0),fT)
 332         call enerprint(energia(0))
 333         write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,49)
 334         write (iout,*) "ftors",ftors
 335 !el        call briefout(i,energia(0))
 336         temp=1.0d0/(beta_h(ib,ipar)*1.987D-3)
 337         write (iout,*) "temp", temp
 338         call pdboutW(i,temp,energia(0),energia(0),0.0d0,0.0d0)
 339 #endif
 340         if (energia(0).ge.1.0d20) then
 341           write (iout,*) "NaNs detected in some of the energy",&
 342            " components for conformation",ii+1
 343           write (iout,*) "The Cartesian geometry is:"
 344           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 345           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 346           write (iout,*) "The internal geometry is:"
 347 !          call intout
 348 !        call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
 349           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 350           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 351           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 352           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 353           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 354           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 355           write (iout,*) "The components of the energy are:"
 356 !          call enerprint(energia(0),fT)
 357           call enerprint(energia(0))
 358           write (iout,*) &
 359             "This conformation WILL NOT be added to the database."
 360           call flush(iout)
 361           goto 121
 362         else
 363 #ifdef DEBUG
 364           if (ipar.eq.iparm) write (iout,*) i,iparm,&
 365             1.0d0/(beta_h(ib,ipar)*1.987D-3),eini,energia(0)
 366 #endif
 367           if (ipar.eq.iparm .and. einicheck.gt.0 .and. &
 368             dabs(eini-energia(0)).gt.tole) then
 369             if (errmsg_count.le.maxerrmsg_count) then
 370               write (iout,'(2a,2e15.5,a,2i8,a,f8.1)') &
 371                "Warning: energy differs remarkably from ",&
 372                " the value read in: ",energia(0),eini," point",&
 373                iii+1,indstart(me1)+iii," T",&
 374                1.0d0/(1.987D-3*beta_h(ib,ipar))
 375 !              call intout
 376               call pdboutW(indstart(me1)+iii,&
 377        1.0d0/(1.987D-3*beta_h(ib,ipar)),&
 378        energia(0),eini,0.0d0,0.0d0)
 379           write (iout,*) "The Cartesian geometry is:"
 380           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 381           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 382           write (iout,*) "The internal geometry is:"
 383 !          call intout
 384 !        call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
 385           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 386           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 387           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 388           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 389           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 390           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 391               call enerprint(energia(0))
 392 !              call enerprint(energia(0),fT)
 393               errmsg_count=errmsg_count+1
 394               if (errmsg_count.gt.maxerrmsg_count) &
 395                 write (iout,*) "Too many warning messages"
 396               if (einicheck.gt.1) then
 397                 write (iout,*) "Calculation stopped."
 398                 call flush(iout)
 399 #ifdef MPI
 400                 call MPI_Abort(WHAM_COMM,IERROR,ERRCODE)
 401 #endif
 402                 call flush(iout)
 403                 return 1
 404               endif
 405             endif
 406           endif
 407           potE(iii+1,iparm)=energia(0)
 408           do k=1,49
 409             enetb(k,iii+1,iparm)=energia(k)
 410           enddo
 411 !           write (iout,'(2i5,21f8.2)') "debug",k,iii+1,(enetb(k,iii+1,iparm),k=1,21)
 412 !           write (iout,*) "debug",k,iii+1,(enetb(k,iii+1,iparm),k=1,21)
 413 #ifdef DEBUG
 414           write (iout,'(2i5,f10.1,3e15.5)') i,iii,&
 415            1.0d0/(beta_h(ib,ipar)*1.987D-3),energia(0),eini,efree
 416 !          call enerprint(energia(0),fT)
 417 #endif
 418 #ifdef DEBUG
 419           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 420           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 421           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 422           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 423           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 424           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 425           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 426           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 427           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 428           write (iout,'(8f10.5)') (q(k,iii+1),k=1,nQ)
 429           write (iout,'(f10.5,i10)') rmsdev,iscor
 430 !          call enerprint(energia(0),fT)
 431           call enerprint(energia(0))
 432         call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
 433 #endif
 434         endif
 435
 436         enddo ! iparm
 437
 438         iii=iii+1
 439         if (q(1,iii).le.0.0d0 .and. indpdb.gt.0) q(1,iii)=qwolynes(0,0)
 440         write (ientout,rec=iii) &
 441          ((csingle(l,k),l=1,3),k=1,nres),&
 442          ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
 443          nss,(ihpb(k),jhpb(k),k=1,nss),&
 444          potE(iii,ipar),efree,rmsdev,(q(k,iii),k=1,nQ),iR,ib,ipar
 445 !        write (iout,'(2i5,2e15.5)') ii,iii,potE(iii,ipar),efree
 446 #ifdef MPI
 447         if (separate_parset) then
 448           snk_p(iR,ib,1)=snk_p(iR,ib,1)+1
 449         else
 450           snk_p(iR,ib,ipar)=snk_p(iR,ib,ipar)+1
 451         endif
 452 !        write (iout,*) "iii",iii," iR",iR," ib",ib," ipar",ipar,
 453 !     &   " snk",snk_p(iR,ib,ipar)
 454 #else
 455         snk(iR,ib,ipar,islice)=snk(iR,ib,ipar,islice)+1
 456 #endif
 457   121   continue
 458       enddo
 459 #ifdef MPI
 460       scount(me)=iii
 461       write (iout,*) "Me",me," scount",scount(me)
 462       call flush(iout)
 463 !  Master gathers updated numbers of conformations written by all procs.
 464       call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount_(0), 1, &
 465         MPI_INTEGER, WHAM_COMM, IERROR)
 466       indstart(0)=1
 467       indend(0)=scount_(0)
 468       do i=1, Nprocs-1
 469         indstart(i)=indend(i-1)+1
 470         indend(i)=indstart(i)+scount_(i)-1
 471       enddo
 472       write (iout,*)
 473       write (iout,*) "Revised conformation counts"
 474       do i=0,nprocs1-1
 475         write (iout,'(a,i5,a,i7,a,i7,a,i7)') &
 476           "Processor",i," indstart",indstart(i),&
 477           " indend",indend(i)," count",scount_(i)
 478       enddo
 479       call flush(iout)
 480       call MPI_AllReduce(snk_p(1,1,1),snk(1,1,1,islice),&
 481         MaxR*MaxT_h*nParmSet,&
 482         MPI_INTEGER,MPI_SUM,WHAM_COMM,IERROR)
 483 #endif
 484       stot(islice)=0
 485       do iparm=1,nParmSet
 486         do ib=1,nT_h(iparm)
 487           do i=1,nR(ib,iparm)
 488             stot(islice)=stot(islice)+snk(i,ib,iparm,islice)
 489           enddo
 490         enddo
 491       enddo
 492       write (iout,*) "Revised SNK"
 493       do iparm=1,nParmSet
 494         do ib=1,nT_h(iparm)
 495           write (iout,'("Param",i3," Temp",f6.1,3x,32i8)') &
 496            iparm,1.0d0/(1.987D-3*beta_h(ib,iparm)),&
 497            (snk(i,ib,iparm,islice),i=1,nR(ib,iparm))
 498           write (iout,*) "snk_p",(snk_p(i,ib,iparm),i=1,nR(ib,iparm))
 499         enddo
 500       enddo
 501       write (iout,'("Total",i10)') stot(islice)
 502       call flush(iout)
 503       do i=0,nprocs
 504         scount(i)=scount_(i)
 505       enddo
 506       return
 507   101 write (iout,*) "Error in scratchfile."
 508       call flush(iout)
 509 !el#undef DEBUG
 510       return 1
 511       end subroutine enecalc
 512 !------------------------------------------------------------------------------
 513       logical function conf_check(ii,iprint)
 514
 515       use geometry_data, only:c,phi,theta,alph,omeg,deg2rad,rad2deg,vbld
 516       use geometry, only:int_from_cart1
 517 !      include "DIMENSIONS"
 518 !      include "DIMENSIONS.ZSCOPT"
 519 !      include "DIMENSIONS.FREE"
 520 #ifdef MPI
 521 !      use MPI_data
 522       include "mpif.h"
 523 !      include "COMMON.MPI"
 524 #endif
 525 !      include "COMMON.CHAIN"
 526 !      include "COMMON.IOUNITS"
 527 !      include "COMMON.PROTFILES"
 528 !      include "COMMON.NAMES"
 529 !      include "COMMON.VAR"
 530 !      include "COMMON.SBRIDGE"
 531 !      include "COMMON.GEO"
 532 !      include "COMMON.FFIELD"
 533 !      include "COMMON.ENEPS"
 534 !      include "COMMON.LOCAL"
 535 !      include "COMMON.WEIGHTS"
 536 !      include "COMMON.INTERACT"
 537 !      include "COMMON.FREE"
 538 !      include "COMMON.ENERGIES"
 539 !      include "COMMON.CONTROL"
 540 !      include "COMMON.TORCNSTR"
 541 !      implicit none
 542 #ifdef MPI
 543       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 544 #endif
 545       integer :: j,k,l,ii,itj,iprint,mnum,mnum1
 546       if (.not.check_conf) then
 547         conf_check=.true.
 548         return
 549       endif
 550       call int_from_cart1(.false.)
 551       do j=nnt+1,nct
 552         mnum=molnum(j)
 553         mnum1=molnum(j-1)
 554         if (mnum.ne.1) cycle
 555         if (itype(j-1,mnum1).ne.ntyp1_molec(mnum1) .and. itype(j,mnum).ne.ntyp1_molec(mnum) .and. &
 556           (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
 557           if (iprint.gt.0) &
 558           write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),&
 559             " for conformation",ii,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
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
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       do iparm=1,iparm
1276         if (iparm.ne.iparmprint) exit
1277         call restore_parm(iparm)
1278         do ib=1,nT_h(iparm)
1279 #ifdef DEBUG
1280           write (iout,*) "iparm",iparm," ib",ib
1281 #endif
1282           temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1283 !          quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1284 !          quotl=1.0d0
1285 !          kfacl=1.0d0
1286 !          do l=1,5
1287 !            quotl1=quotl
1288 !            quotl=quotl*quot
1289 !            kfacl=kfacl*kfac
1290 !            fT(l)=kfacl/(kfacl-1.0d0+quotl)
1291 !          enddo
1292 !el old rescale weights
1293 !
1294 !            if (rescale_mode.eq.1) then
1295 !              quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1296 #if defined(FUNCTH)
1297               tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1298               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
1299 #elif defined(FUNCT)
1300               ft(6)=quot
1301 #else
1302               ft(6)=1.0d0
1303 #endif
1304 !              quotl=1.0d0
1305 !              kfacl=1.0d0
1306 !              do l=1,5
1307 !                quotl1=quotl
1308 !                quotl=quotl*quot
1309 !                kfacl=kfacl*kfac
1310 !                fT(l)=kfacl/(kfacl-1.0d0+quotl)
1311 !              enddo
1312 !            else if (rescale_mode.eq.2) then
1313 !              quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1314 !#if defined(FUNCTH)
1315 !              tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1316 !              ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0
1317 !#elif defined(FUNCT)
1318 !              ft(6)=quot
1319 !#else
1320 !              ft(6)=1.0d0
1321 !#endif
1322 !              quotl=1.0d0
1323 !              do l=1,5
1324 !                quotl=quotl*quot
1325 !                fT(l)=1.12692801104297249644d0/ &
1326 !                   dlog(dexp(quotl)+dexp(-quotl))
1327 !              enddo
1328 !              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
1329 !            else if (rescale_mode.eq.0) then
1330 !              do l=1,5
1331 !                fT(l)=0.0d0
1332 !              enddo
1333 !            else
1334 !              write (iout,*) &
1335 !              "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
1336 !              call flush(iout)
1337 !              return 1
1338 !            endif
1339 ! el end old rescale weihgts
1340           call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3))
1341
1342 #ifdef MPI
1343           do i=1,scount(me1)
1344 #else
1345           do i=1,ntot(islice)
1346 #endif
1347             evdw=enetb(1,i,iparm)
1348 !            evdw_t=enetb(21,i,iparm)
1349             evdw_t=enetb(20,i,iparm)
1350 #ifdef SCP14
1351 !            evdw2_14=enetb(17,i,iparm)
1352             evdw2_14=enetb(18,i,iparm)
1353             evdw2=enetb(2,i,iparm)+evdw2_14
1354 #else
1355             evdw2=enetb(2,i,iparm)
1356             evdw2_14=0.0d0
1357 #endif
1358 #ifdef SPLITELE
1359             ees=enetb(3,i,iparm)
1360             evdw1=enetb(16,i,iparm)
1361 #else
1362             ees=enetb(3,i,iparm)
1363             evdw1=0.0d0
1364 #endif
1365             ecorr=enetb(4,i,iparm)
1366             ecorr5=enetb(5,i,iparm)
1367             ecorr6=enetb(6,i,iparm)
1368             eel_loc=enetb(7,i,iparm)
1369             eello_turn3=enetb(8,i,iparm)
1370             eello_turn4=enetb(9,i,iparm)
1371             eello_turn6=enetb(10,i,iparm)
1372             ebe=enetb(11,i,iparm)
1373             escloc=enetb(12,i,iparm)
1374             etors=enetb(13,i,iparm)
1375             etors_d=enetb(14,i,iparm)
1376             ehpb=enetb(15,i,iparm)
1377
1378             estr=enetb(17,i,iparm)
1379 !            estr=enetb(18,i,iparm)
1380 !            esccor=enetb(19,i,iparm)
1381             esccor=enetb(21,i,iparm)
1382 !            edihcnstr=enetb(20,i,iparm)
1383             edihcnstr=enetb(19,i,iparm)
1384             ecationcation=enetb(42,i,iparm)
1385             ecation_prot=enetb(41,i,iparm)
1386             evdwpp =      enetb(26,i,iparm)
1387             eespp =      enetb(27,i,iparm)
1388             evdwpsb =      enetb(28,i,iparm)
1389             eelpsb =      enetb(29,i,iparm)
1390             evdwsb =      enetb(30,i,iparm)
1391             eelsb =      enetb(31,i,iparm)
1392             estr_nucl =      enetb(32,i,iparm)
1393             ebe_nucl =      enetb(33,i,iparm)
1394             esbloc =      enetb(34,i,iparm)
1395             etors_nucl =      enetb(35,i,iparm)
1396             etors_d_nucl =      enetb(36,i,iparm)
1397             ecorr_nucl =      enetb(37,i,iparm)
1398             ecorr3_nucl =      enetb(38,i,iparm)
1399             epeppho=   enetb(49,i,iparm)
1400             escpho=    enetb(48,i,iparm)
1401             epepbase=  enetb(47,i,iparm)
1402             escbase=   enetb(46,i,iparm)
1403 !      wscbase=ww(46)
1404 !      wpepbase=ww(47)
1405 !      wscpho=ww(48)
1406 !      wpeppho=ww(49)
1407
1408 !#ifdef SPLITELE
1409 !            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1410 !            +wvdwpp*evdw1 &
1411 !            +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1412 !            +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1413 !            +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1414 !            +ft(2)*wturn3*eello_turn3 &
1415 !            +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc &
1416 !            +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1417 !            +wbond*estr
1418 !#else
1419 !            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1420 !            +ft(1)*welec*(ees+evdw1) &
1421 !            +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1422 !            +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1423 !            +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1424 !            +ft(2)*wturn3*eello_turn3 &
1425 !            +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1426 !            +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1427 !            +wbond*estr
1428 !#endif
1429
1430 #ifdef SPLITELE
1431             etot=wsc*evdw+wscp*evdw2+welec*ees &
1432             +wvdwpp*evdw1 &
1433             +wang*ebe+wtor*etors+wscloc*escloc &
1434             +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
1435             +wcorr6*ecorr6+wturn4*eello_turn4 &
1436             +wturn3*eello_turn3 &
1437             +wturn6*eello_turn6+wel_loc*eel_loc &
1438             +edihcnstr+wtor_d*etors_d+wsccor*esccor &
1439             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
1440             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1441             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1442             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1443             +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
1444             +wscbase*escbase&
1445             +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
1446
1447 #else
1448             etot=wsc*evdw+wscp*evdw2 &
1449             +welec*(ees+evdw1) &
1450             +wang*ebe+wtor*etors+wscloc*escloc &
1451             +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
1452             +wcorr6*ecorr6+wturn4*eello_turn4 &
1453             +wturn3*eello_turn3 &
1454             +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
1455             +wtor_d*etors_d+wsccor*esccor &
1456             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
1457             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1458             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1459             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1460             +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
1461             +wscbase*escbase&
1462             +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
1463
1464 #endif
1465
1466 #ifdef MPI
1467             Fdimless(i)= &
1468               beta_h(ib,iparm)*etot-entfac(i)
1469             potE(i,iparm)=etot
1470 #ifdef DEBUG
1471             write (iout,*) i,indstart(me)+i-1,ib,&
1472              1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),&
1473              -entfac(i),Fdimless(i)
1474 #endif
1475 #else
1476             Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
1477             potE(i,iparm)=etot
1478 #endif
1479           enddo   ! i
1480 #ifdef MPI
1481           scount_(:)=0
1482           do i=1,scount(me1)
1483             Fdimless_(i)=Fdimless(i)
1484           enddo
1485           write(iout,*) "before gather Fdimless"
1486           write(iout,*) scount(me),scount_(0),idispl(0)
1487           call MPI_Gatherv(Fdimless_(1),scount(me),&
1488            MPI_REAL,Fdimless(1),&
1489            scount_(0),idispl(0),MPI_REAL,Master,&
1490            WHAM_COMM, IERROR)
1491           write(iout,*) "after gather Fdimless"
1492 #ifdef DEBUG
1493           call MPI_Gatherv(potE(1,iparm),scount_(me),&
1494            MPI_DOUBLE_PRECISION,potE(1,iparm),&
1495            scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1496            WHAM_COMM, IERROR)
1497           call MPI_Gatherv(entfac(1),scount(me),&
1498            MPI_DOUBLE_PRECISION,entfac(1),&
1499            scount_(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1500            WHAM_COMM, IERROR)
1501 #endif
1502           if (me.eq.Master) then
1503 #ifdef DEBUG
1504           write (iout,*) "The FDIMLESS array before sorting"
1505           do i=1,ntot(islice)
1506             write (iout,*) i,fdimless(i)
1507           enddo
1508 #endif
1509 #endif
1510           do i=1,ntot(islice)
1511             iperm(i)=i
1512           enddo
1513           call mysort1(ntot(islice),Fdimless,iperm)
1514 #ifdef DEBUG
1515           write (iout,*) "The FDIMLESS array after sorting"
1516           do i=1,ntot(islice)
1517             write (iout,*) i,iperm(i),fdimless(i)
1518           enddo
1519 #endif
1520           qfree=0.0d0
1521           do i=1,ntot(islice)
1522             qfree=qfree+exp(-fdimless(i)+fdimless(1))
1523           enddo
1524 !          write (iout,*) "qfree",qfree
1525           nlist=1
1526           sumprob=0.0
1527           do i=1,min0(ntot(islice),ensembles)
1528             sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
1529 #ifdef DEBUG
1530             write (iout,*) i,ib,beta_h(ib,iparm),&
1531              1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),&
1532              potE(iperm(i),iparm),&
1533              -entfac(iperm(i)),fdimless(i),sumprob
1534 #endif
1535             if (sumprob.gt.0.99d0) goto 122
1536             nlist=nlist+1
1537           enddo
1538   122     continue
1539 #ifdef MPI
1540           endif
1541           call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,&
1542              IERROR)
1543           call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,&
1544              IERROR)
1545           do i=1,nlist
1546             ii=iperm(i)
1547             iproc=0
1548             do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
1549               iproc=iproc+1
1550             enddo
1551             if (iproc.ge.nprocs) then
1552               write (iout,*) "Fatal error: processor out of range",iproc
1553               call flush(iout)
1554               if (bxfile) then
1555                 close (ientout)
1556               else
1557                 close (ientout,status="delete")
1558               endif
1559               return 1
1560             endif
1561             ik=ii-indstart(iproc)+1
1562             if (iproc.ne.Master) then
1563               if (me.eq.iproc) then
1564 #ifdef DEBUG
1565                 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,&
1566                  " energy",potE(ik,iparm)
1567 #endif
1568                 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,&
1569                   Master,i,WHAM_COMM,IERROR)
1570               else if (me.eq.Master) then
1571                 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,&
1572                   WHAM_COMM,STATUS,IERROR)
1573               endif
1574             else if (me.eq.Master) then
1575               enepot(i)=potE(ik,iparm)
1576             endif
1577           enddo
1578 #else
1579           do i=1,nlist
1580             enepot(i)=potE(iperm(i),iparm)
1581           enddo
1582 #endif
1583 #ifdef MPI
1584           if (me.eq.Master) then
1585 #endif
1586           write(licz3,'(bz,i3.3)') iparm
1587           write(licz2,'(bz,i2.2)') islice
1588           if (temper.lt.100.0d0) then
1589             write(ctemper,'(f3.0)') temper
1590           else if (temper.lt.1000.0) then
1591             write (ctemper,'(f4.0)') temper
1592           else
1593             write (ctemper,'(f5.0)') temper
1594           endif
1595           if (nparmset.eq.1) then
1596             if (separate_parset) then
1597               write(licz4,'(bz,i3.3)') myparm
1598               pdbname=prefix(:ilen(prefix))//"_par"//licz4
1599             else
1600               pdbname=prefix(:ilen(prefix))
1601             endif
1602           else
1603             pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
1604           endif
1605           if (nslice.eq.1) then
1606             pdbname=pdbname(:ilen(pdbname))//"_T_"// &
1607               ctemper(:ilen(ctemper))//"pdb"
1608           else
1609             pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"// &
1610               ctemper(:ilen(ctemper))//"pdb"
1611           endif
1612           open(ipdb,file=pdbname)
1613           do i=1,nlist
1614             read (ientout,rec=iperm(i)) &
1615               ((csingle(l,k),l=1,3),k=1,nres),&
1616               ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1617               nss,(ihpb(k),jhpb(k),k=1,nss),&
1618               eini,efree,rmsdev,iscor
1619             do j=1,2*nres
1620               do k=1,3
1621                 c(k,j)=csingle(k,j)
1622               enddo
1623             enddo
1624             call returnbox
1625             do j=1,2*nres
1626               do k=1,3
1627                 csingle(k,j)=c(k,j)
1628               enddo
1629             enddo
1630
1631             eini=fdimless(i)
1632             call pdboutW(iperm(i),temper,eini,enepot(i),efree,rmsdev)
1633           enddo
1634 #ifdef MPI
1635           endif
1636 #endif
1637         enddo     ! ib
1638       enddo       ! iparm
1639       if (bxfile) then
1640         close(ientout)
1641       else
1642         close(ientout,status="delete")
1643       endif
1644       do i=0,nprocs
1645         scount(i)=scount_(i)
1646       enddo
1647       return
1648       end subroutine make_ensembles
1649 !--------------------------------------------------------------------------
1650       subroutine mysort1(n, x, ipermut)
1651 !      implicit none
1652       integer :: i,j,imax,ipm,n
1653       real(kind=4) :: x(n)
1654       integer :: ipermut(n)
1655       real(kind=4) :: xtemp
1656       do i=1,n
1657         xtemp=x(i)
1658         imax=i
1659         do j=i+1,n
1660           if (x(j).lt.xtemp) then
1661             imax=j
1662             xtemp=x(j)
1663           endif
1664         enddo
1665         x(imax)=x(i)
1666         x(i)=xtemp
1667         ipm=ipermut(imax)
1668         ipermut(imax)=ipermut(i)
1669         ipermut(i)=ipm
1670       enddo
1671       return
1672       end subroutine mysort1
1673 !--------------------------------------------------------------------------
1674       subroutine alloc_enecalc_arrays(iparm)
1675
1676       use control_data
1677       use geometry_data, only:maxlob
1678       integer :: iparm
1679 !---------------------------
1680 ! COMMON.ENERGIES form wham_data
1681 !      common /energies/
1682       allocate(potE(MaxStr_Proc,iparm)) !(MaxStr_Proc,Max_Parm)
1683       allocate(entfac(MaxStr_Proc)) !(MaxStr_Proc)
1684       allocate(q(nQ+2,MaxStr_Proc)) !(MaxQ+2,MaxStr_Proc)
1685       allocate(enetb(max_ene,MaxStr_Proc,iparm)) !(max_ene,MaxStr_Proc,Max_Parm)
1686 !
1687 ! allocate ENECALC arrays
1688 !---------------------------
1689 ! COMMON.ALLPARM
1690 !      common /allparm/
1691       allocate(ww_all(max_eneW,iparm)) !(max_ene,max_parm) ! max_eneW
1692       allocate(vbldp0_all(iparm),akp_all(nParmSet)) !(max_parm)
1693       allocate(vbldsc0_all(maxbondterm,ntyp,iparm),&
1694         aksc_all(maxbondterm,ntyp,iparm),&
1695         abond0_all(maxbondterm,ntyp,iparm)) !(maxbondterm,ntyp,max_parm)
1696       allocate(a0thet_all(-ntyp:ntyp,iparm)) !(-ntyp:ntyp,max_parm)
1697       allocate(athet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm),&
1698         bthet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm)) !(2,-ntyp:ntyp,-1:1,-1:1,max_parm)
1699       allocate(polthet_all(0:3,-ntyp:ntyp,iparm)) !(0:3,-ntyp:ntyp,max_parm)
1700       allocate(gthet_all(3,-ntyp:ntyp,iparm)) !(3,-ntyp:ntyp,max_parm)
1701       allocate(theta0_all(-ntyp:ntyp,iparm),&
1702         sig0_all(-ntyp:ntyp,iparm),sigc0_all(-ntyp:ntyp,nParmSet)) !(-ntyp:ntyp,max_parm)
1703       allocate(aa0thet_all(-nthetyp-1:nthetyp+1,&
1704         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1705 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1706       allocate(aathet_all(ntheterm,-nthetyp-1:nthetyp+1,&
1707         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1708 !(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1709       allocate(bbthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1710         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1711       allocate(ccthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1712         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1713       allocate(ddthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1714         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1715       allocate(eethet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1716         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1717 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,
1718 !     & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1719       allocate(ffthet_all1(ndouble,ndouble,ntheterm3,&
1720         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1721         -nthetyp-1:nthetyp+1,iparm))
1722       allocate(ggthet_all1(ndouble,ndouble,ntheterm3,&
1723         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1724         -nthetyp-1:nthetyp+1,iparm))
1725       allocate(ffthet_all2(ndouble,ndouble,ntheterm3,&
1726         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1727         -nthetyp-1:nthetyp+1,iparm))
1728       allocate(ggthet_all2(ndouble,ndouble,ntheterm3,&
1729         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1730         -nthetyp-1:nthetyp+1,iparm))
1731 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1732 !-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm)
1733       allocate(dsc_all(ntyp1,iparm),dsc0_all(ntyp1,nParmSet)) !(ntyp1,max_parm)
1734       allocate(bsc_all(maxlob,ntyp,iparm))
1735 !(maxlob,ntyp,max_parm)
1736       allocate(censc_all(3,maxlob,-ntyp:ntyp,iparm)) !(3,maxlob,-ntyp:ntyp,max_parm)
1737       allocate(gaussc_all(3,3,maxlob,-ntyp:ntyp,iparm)) !(3,3,maxlob,-ntyp:ntyp,max_parm)
1738       allocate(sc_parmin_all(65,ntyp,iparm)) !(65,ntyp,max_parm)
1739       allocate(v0_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1740 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1741       allocate(v1_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1742       allocate(v2_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1743 !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1744       allocate(vlor1_all(maxlor,ntortyp,ntortyp,iparm))
1745       allocate(vlor2_all(maxlor,ntortyp,ntortyp,iparm))
1746       allocate(vlor3_all(maxlor,ntortyp,ntortyp,iparm)) !(maxlor,maxtor,maxtor,max_parm)
1747       allocate(v1c_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1748         -ntortyp:ntortyp,2,iparm))
1749       allocate(v1s_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1750         -ntortyp:ntortyp,2,iparm))
1751 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1752       allocate(v2c_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1753         -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1754 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1755       allocate(v2s_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1756         -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1757 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1758       allocate(b1_all(2,-nloctyp:nloctyp,iparm))
1759       allocate(b2_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1760       allocate(cc_all(2,2,-nloctyp:nloctyp,iparm))
1761       allocate(dd_all(2,2,-nloctyp:nloctyp,iparm))
1762       allocate(ee_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1763       allocate(ctilde_all(2,2,-nloctyp:nloctyp,iparm))
1764       allocate(dtilde_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1765       allocate(b1tilde_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1766       allocate(app_all(2,2,iparm),bpp_all(2,2,nParmSet),&
1767         ael6_all(2,2,iparm),ael3_all(2,2,nParmSet)) !(2,2,max_parm)
1768       allocate(aad_all(ntyp,2,iparm),bad_all(ntyp,2,nParmSet)) !(ntyp,2,max_parm)
1769       allocate(aa_all(ntyp,ntyp,iparm),bb_all(ntyp,ntyp,nParmSet),&
1770         augm_all(ntyp,ntyp,iparm),eps_all(ntyp,ntyp,nParmSet),&
1771         sigma_all(ntyp,ntyp,iparm),r0_all(ntyp,ntyp,nParmSet),&
1772         chi_all(ntyp,ntyp,iparm)) !(ntyp,ntyp,max_parm)
1773       allocate(chip_all(ntyp,iparm),alp_all(ntyp,nParmSet)) !(ntyp,max_parm)
1774       allocate(ebr_all(iparm),d0cm_all(nParmSet),akcm_all(nParmSet),&
1775         akth_all(iparm),akct_all(nParmSet),v1ss_all(nParmSet),&
1776         v2ss_all(iparm),v3ss_all(nParmSet)) !(max_parm)
1777       allocate(v1sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1778       allocate(v2sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1779 !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
1780       allocate(nlob_all(ntyp1,iparm)) !(ntyp1,max_parm)
1781       allocate(nlor_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1782       allocate(nterm_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1783 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1784       allocate(ntermd1_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1785         -ntortyp:ntortyp,2,iparm))
1786       allocate(ntermd2_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1787         -ntortyp:ntortyp,2,iparm))
1788 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1789       allocate(nbondterm_all(ntyp,iparm)) !(ntyp,max_parm)
1790       allocate(ithetyp_all(-ntyp1:ntyp1,iparm)) !(-ntyp1:ntyp1,max_parm)
1791       allocate(nthetyp_all(iparm),ntheterm_all(nParmSet),&
1792         ntheterm2_all(iparm),ntheterm3_all(nParmSet),&
1793         nsingle_all(iparm),&
1794         ndouble_all(iparm),nntheterm_all(nParmSet)) !(max_parm)
1795       allocate(nterm_sccor_all(-ntyp:ntyp,-ntyp:ntyp,iparm)) !(-ntyp:ntyp,-ntyp:ntyp,max_parm)
1796 !
1797       end subroutine alloc_enecalc_arrays
1798 !--------------------------------------------------------------------------
1799 !--------------------------------------------------------------------------
1800       end module ene_calc