source/wham/src-M-SAXS.safe/wham_calc1.F

   1       subroutine WHAM_CALC(islice,*)
   2 ! Weighed Histogram Analysis Method (WHAM) code
   3 ! Written by A. Liwo based on the work of Kumar et al.,
   4 ! J.Comput.Chem., 13, 1011 (1992)
   5 !
   6 ! 2/1/05 Multiple temperatures allowed.
   7 ! 2/2/05 Free energies calculated directly from data points
   8 !  acc. to Eq. (21) of Kumar et al.; final histograms also
   9 !  constructed based on this equation.
  10 ! 2/12/05 Multiple parameter sets included
  11 !
  12 ! 2/2/05 Parallel version
  13       implicit none
  14       include "DIMENSIONS"
  15       include "DIMENSIONS.ZSCOPT"
  16       include "DIMENSIONS.FREE"
  17       integer nGridT
  18       parameter (NGridT=400)
  19       integer MaxBinRms,MaxBinRgy
  20       parameter (MaxBinRms=100,MaxBinRgy=100)
  21 c      integer MaxHdim
  22 c      parameter (MaxHdim=200)
  23       integer maxinde
  24       parameter (maxinde=200)
  25 #ifdef MPI
  26       include "mpif.h"
  27       include "COMMON.MPI"
  28       integer ierror,errcode,status(MPI_STATUS_SIZE)
  29 #endif
  30       include "COMMON.CONTROL"
  31       include "COMMON.IOUNITS"
  32       include "COMMON.FREE"
  33       include "COMMON.ENERGIES"
  34       include "COMMON.FFIELD"
  35       include "COMMON.SBRIDGE"
  36       include "COMMON.PROT"
  37       include "COMMON.ENEPS"
  38       include "COMMON.SHIELD"
  39       integer MaxPoint,MaxPointProc
  40       parameter (MaxPoint=MaxStr,
  41      & MaxPointProc=MaxStr_Proc)
  42       double precision finorm_max,potfac,entmin,entmax,expfac,vf
  43       parameter (finorm_max=1.0d0)
  44       integer islice
  45       integer i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln
  46       integer start,end,iharm,ib,iib,nbin1,nbin,nbin_rms,nbin_rgy,
  47      &  nbin_rmsrgy,liczba,iparm,nFi,indrgy,indrms
  48       integer htot(0:MaxHdim),histent(0:2000)
  49       double precision v(MaxPointProc,MaxR,MaxT_h,Max_Parm)
  50       double precision energia(0:max_ene)
  51 #ifdef MPI
  52       integer tmax_t,upindE_p
  53       double precision fi_p(MaxR,MaxT_h,Max_Parm),
  54      & fimax_p(MaxR,MaxT_h,Max_Parm)
  55       double precision sumW_p(0:nGridT,Max_Parm),
  56      & sumE_p(0:nGridT,Max_Parm),sumEsq_p(0:nGridT,Max_Parm),
  57      & sumQ_p(MaxQ1,0:nGridT,Max_Parm),
  58      & sumQsq_p(MaxQ1,0:nGridT,Max_Parm),
  59      & sumEQ_p(MaxQ1,0:nGridT,Max_Parm),
  60      & sumEprim_p(MaxQ1,0:nGridT,Max_Parm),
  61      & sumEbis_p(0:nGridT,Max_Parm)
  62       double precision hfin_p(0:MaxHdim,maxT_h),
  63      & hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,
  64      & hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
  65       double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
  66       double precision potEmin_t,entmin_p,entmax_p
  67       double precision ePMF,ePMF_q
  68       double precision weimax_(0:ngridT)
  69       integer histent_p(0:2000)
  70       logical lprint /.true./
  71 #endif
  72       double precision delta_T /1.0d0/
  73       double precision rgymin,rmsmin,rgymax,rmsmax
  74       double precision sumW(0:NGridT,Max_Parm),sumE(0:NGridT,Max_Parm),
  75      & sumEsq(0:NGridT,Max_Parm),sumQ(MaxQ1,0:NGridT,Max_Parm),
  76      & sumQsq(MaxQ1,0:NGridT,Max_Parm),sumEQ(MaxQ1,0:NGridT,Max_Parm),
  77      & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT,
  78      & weight,econstr
  79       double precision fi(MaxR,maxT_h,Max_Parm),
  80      & fimax(MaxR,maxT_h,Max_Parm),
  81      & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,
  82      & hfin(0:MaxHdim,maxT_h),histE(0:maxindE),
  83      & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),
  84      & potEmin,ent,
  85      & hfin_ent(0:MaxHdim),vmax,aux,weimax(0:nGridT,Max_Parm)
  86       double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
  87      &  eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
  88      &  eplus,eminus,logfac,tanhT,tt
  89       double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
  90      &  escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
  91      &  eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
  92      &  eliptran,esaxs
  93
  94       integer ind_point(maxpoint),upindE,indE
  95       character*16 plik
  96       character*1 licz1
  97       character*2 licz2
  98       character*3 licz3
  99       character*128 nazwa
 100       integer ilen
 101       external ilen
 102
 103       write(licz2,'(bz,i2.2)') islice
 104       nbin1 = 1.0d0/delta
 105       write (iout,'(//80(1h-)/"Solving WHAM equations for slice",
 106      &  i2/80(1h-)//)') islice
 107       write (iout,*) "delta",delta," nbin1",nbin1
 108       write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
 109       call flush(iout)
 110       dmin=0.0d0
 111       tmax=0
 112       potEmin=1.0d10
 113       rgymin=1.0d10
 114       rmsmin=1.0d10
 115       rgymax=0.0d0
 116       rmsmax=0.0d0
 117       do t=0,MaxN
 118         htot(t)=0
 119       enddo
 120 C#define DEBUG
 121 #ifdef MPI
 122       do i=1,scount(me1)
 123 #else
 124       do i=1,ntot(islice)
 125 #endif
 126 c        write (iout,*) "i",i," potE",(potE(i,j),j=1,nParmset)
 127         do j=1,nParmSet
 128           if (potE(i,j).le.potEmin) potEmin=potE(i,j)
 129         enddo
 130         if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
 131         if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
 132         if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
 133         if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
 134         ind_point(i)=0
 135         do j=nQ,1,-1
 136           ind=(q(j,i)-dmin+1.0d-8)/delta
 137           if (j.eq.1) then
 138             ind_point(i)=ind_point(i)+ind
 139           else
 140             ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
 141           endif
 142 c          write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
 143 c     &      ind_point(i)
 144 c          call flush(iout)
 145           if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
 146             write (iout,*) "Error - index exceeds range for point",i,
 147      &      " q=",q(j,i)," ind",ind_point(i)
 148 #ifdef MPI
 149             write (iout,*) "Processor",me1
 150             call flush(iout)
 151             call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
 152 #endif
 153             stop
 154           endif
 155         enddo ! j
 156         if (ind_point(i).gt.tmax) tmax=ind_point(i)
 157         htot(ind_point(i))=htot(ind_point(i))+1
 158 #ifdef DEBUG
 159         write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),
 160      &   " htot",htot(ind_point(i))
 161         call flush(iout)
 162 #endif
 163       enddo ! i
 164
 165       write (iout,*) "potEmin before reduce",potEmin
 166       nbin=nbin1**nQ-1
 167       write (iout,'(a)') "Numbers of counts in Q bins"
 168       do t=0,tmax
 169         if (htot(t).gt.0) then
 170         write (iout,'(i15,$)') t
 171         liczba=t
 172         do j=1,nQ
 173           jj = mod(liczba,nbin1)
 174           liczba=liczba/nbin1
 175           write (iout,'(i5,$)') jj
 176         enddo
 177         write (iout,'(i8)') htot(t)
 178         endif
 179       enddo
 180       do iparm=1,nParmSet
 181       write (iout,'(a,i3)') "Number of data points for parameter set",
 182      & iparm
 183       write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),
 184      & ib=1,nT_h(iparm))
 185       write (iout,'(i8)') stot(islice)
 186       write (iout,'(a)')
 187       enddo
 188       call flush(iout)
 189
 190 #ifdef MPI
 191       call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
 192      &  WHAM_COMM,IERROR)
 193       tmax=tmax_t
 194       call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,
 195      &  MPI_MIN,WHAM_COMM,IERROR)
 196       call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,
 197      &  MPI_MIN,WHAM_COMM,IERROR)
 198       call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,
 199      &  MPI_MAX,WHAM_COMM,IERROR)
 200       call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,
 201      &  MPI_MIN,WHAM_COMM,IERROR)
 202       call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
 203      &  MPI_MAX,WHAM_COMM,IERROR)
 204 c      potEmin=potEmin_t/2
 205       potEmin=potEmin_t
 206       rgymin=rgymin_t
 207       rgymax=rgymax_t
 208       rmsmin=rmsmin_t
 209       rmsmax=rmsmax_t
 210       write (iout,*) "potEmin",potEmin
 211 #endif
 212       rmsmin=deltrms*dint(rmsmin/deltrms)
 213       rmsmax=deltrms*dint(rmsmax/deltrms)
 214       rgymin=deltrms*dint(rgymin/deltrgy)
 215       rgymax=deltrms*dint(rgymax/deltrgy)
 216       nbin_rms=(rmsmax-rmsmin)/deltrms
 217       nbin_rgy=(rgymax-rgymin)/deltrgy
 218       write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,
 219      & " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
 220       nFi=0
 221       do i=1,nParmSet
 222         do j=1,nT_h(i)
 223           nFi=nFi+nR(j,i)
 224         enddo
 225       enddo
 226       write (iout,*) "nFi",nFi
 227 ! Compute the Boltzmann factor corresponing to restrain potentials in different
 228 ! simulations.
 229 #ifdef MPI
 230       do i=1,scount(me1)
 231 #else
 232       do i=1,ntot(islice)
 233 #endif
 234 c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
 235         do iparm=1,nParmSet
 236 #ifdef DEBUG
 237           write (iout,'(2i5,21f8.2)') i,iparm,
 238      &     (enetb(k,i,iparm),k=1,22)
 239 #endif
 240           call restore_parm(iparm)
 241 #ifdef DEBUG
 242           write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
 243      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
 244      &      wtor_d,wsccor,wbond
 245 #endif
 246           do ib=1,nT_h(iparm)
 247             if (rescale_mode.eq.1) then
 248               quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 249               quotl=1.0d0
 250               kfacl=1.0d0
 251               do l=1,5
 252                 quotl1=quotl
 253                 quotl=quotl*quot
 254                 kfacl=kfacl*kfac
 255                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
 256               enddo
 257 #if defined(FUNCTH)
 258               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 259               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 260 #elif defined(FUNCT)
 261               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 262 #else
 263               ft(6)=1.0d0
 264 #endif
 265             else if (rescale_mode.eq.2) then
 266               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 267               quotl=1.0d0
 268               do l=1,5
 269                 quotl=quotl*quot
 270                 fT(l)=1.12692801104297249644d0/
 271      &             dlog(dexp(quotl)+dexp(-quotl))
 272               enddo
 273 #if defined(FUNCTH)
 274               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 275               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 276 #elif defined(FUNCT)
 277               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 278 #else
 279               ft(6)=1.0d0
 280 #endif
 281 c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 282             else if (rescale_mode.eq.0) then
 283               do l=1,6
 284                 fT(l)=1.0d0
 285               enddo
 286             else
 287               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
 288      &          rescale_mode
 289               call flush(iout)
 290               return1
 291             endif
 292             evdw=enetb(1,i,iparm)
 293             evdw_t=enetb(21,i,iparm)
 294 #ifdef SCP14
 295             evdw2_14=enetb(17,i,iparm)
 296             evdw2=enetb(2,i,iparm)+evdw2_14
 297 #else
 298             evdw2=enetb(2,i,iparm)
 299             evdw2_14=0.0d0
 300 #endif
 301 #ifdef SPLITELE
 302             ees=enetb(3,i,iparm)
 303             evdw1=enetb(16,i,iparm)
 304 #else
 305             ees=enetb(3,i,iparm)
 306             evdw1=0.0d0
 307 #endif
 308             ecorr=enetb(4,i,iparm)
 309             ecorr5=enetb(5,i,iparm)
 310             ecorr6=enetb(6,i,iparm)
 311             eel_loc=enetb(7,i,iparm)
 312             eello_turn3=enetb(8,i,iparm)
 313             eello_turn4=enetb(9,i,iparm)
 314             eturn6=enetb(10,i,iparm)
 315             ebe=enetb(11,i,iparm)
 316             escloc=enetb(12,i,iparm)
 317             etors=enetb(13,i,iparm)
 318             etors_d=enetb(14,i,iparm)
 319             ehpb=enetb(15,i,iparm)
 320             estr=enetb(18,i,iparm)
 321             esccor=enetb(19,i,iparm)
 322             edihcnstr=enetb(20,i,iparm)
 323             eliptran=enetb(22,i,iparm)
 324             esaxs=enetb(26,i,iparm)
 325
 326 #ifdef DEBUG
 327             write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
 328      &       evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
 329      &       etors,etors_d,eello_turn3,eello_turn4,esccor,esaxs
 330 #endif
 331
 332 #ifdef SPLITELE
 333             if (shield_mode.gt.0) then
 334             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 335      &      +ft(1)*welec*ees
 336      &      +ft(1)*wvdwpp*evdw1
 337      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 338 !     &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 339      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 340      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 341      &      +ft(2)*wturn3*eello_turn3
 342      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 343      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 344      &      +wbond*estr+wliptran*eliptran+wsaxs*esaxs
 345              else
 346             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
 347      &      +wvdwpp*evdw1
 348      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 349 c     &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 350      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 351      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 352      &      +ft(2)*wturn3*eello_turn3
 353      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 354      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 355      &      +wbond*estr+wliptran*eliptran+wsaxs*esaxs
 356              endif
 357 #else
 358       if (shield_mode.gt.0) then
 359             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 360      &      +ft(1)*welec*(ees+evdw1)
 361      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 362 c     &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 363      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 364      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 365      &      +ft(2)*wturn3*eello_turn3
 366      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 367      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 368      &      +wbond*estr+wliptran*eliptran+wsaxs*esaxs
 369             else
 370             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 371      &      +ft(1)*welec*(ees+evdw1)
 372      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 373 c     &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 374      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 375      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 376      &      +ft(2)*wturn3*eello_turn3
 377      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 378      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 379      &      +wbond*estr+wliptran*eliptran+wsaxs*esaxs
 380            endif
 381
 382 #endif
 383 #ifdef DEBUG
 384             write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
 385      &        etot,potEmin
 386 #endif
 387 #ifdef DEBUG
 388             if (iparm.eq.1 .and. ib.eq.1) then
 389             write (iout,*)"Conformation",i
 390             energia(0)=etot
 391             do k=1,max_ene
 392               energia(k)=enetb(k,i,iparm)
 393             enddo
 394             call enerprint(energia(0),fT)
 395             endif
 396 #endif
 397             do kk=1,nR(ib,iparm)
 398               Econstr=0.0d0
 399               do j=1,nQ
 400                 dd = q(j,i)
 401                 Econstr=Econstr+Kh(j,kk,ib,iparm)
 402      &           *(dd-q0(j,kk,ib,iparm))**2
 403               enddo
 404 c Adaptive potential contribution
 405               if (adaptive) then
 406                 call PMF_energy(q(1,i),ib,kk,iparm,ePMF,ePMF_q)
 407                 Econstr=Econstr+ePMF
 408               endif
 409               v(i,kk,ib,iparm)=
 410      &          -beta_h(ib,iparm)*(etot-potEmin+Econstr)
 411 #ifdef DEBUG
 412               write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
 413      &         etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
 414 #endif
 415             enddo ! kk
 416           enddo   ! ib
 417         enddo     ! iparm
 418       enddo       ! i
 419 ! Simple iteration to calculate free energies corresponding to all simulation
 420 ! runs.
 421       do iter=1,maxit
 422
 423 ! Compute new free-energy values corresponding to the righ-hand side of the
 424 ! equation and their derivatives.
 425         write (iout,*) "------------------------fi"
 426 #ifdef MPI
 427         do t=1,scount(me1)
 428 #else
 429         do t=1,ntot(islice)
 430 #endif
 431           vmax=-1.0d+20
 432           do i=1,nParmSet
 433             do k=1,nT_h(i)
 434               do l=1,nR(k,i)
 435                 vf=v(t,l,k,i)+f(l,k,i)
 436                 if (vf.gt.vmax) vmax=vf
 437               enddo
 438             enddo
 439           enddo
 440           denom=0.0d0
 441           do i=1,nParmSet
 442             do k=1,nT_h(i)
 443               do l=1,nR(k,i)
 444                 aux=f(l,k,i)+v(t,l,k,i)-vmax
 445                 if (aux.gt.-200.0d0)
 446      &          denom=denom+snk(l,k,i,islice)*dexp(aux)
 447               enddo
 448             enddo
 449           enddo
 450           entfac(t)=-dlog(denom)-vmax
 451 #ifdef DEBUG
 452           write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
 453 #endif
 454         enddo
 455
 456         do iparm=1,nParmSet
 457           do iib=1,nT_h(iparm)
 458             do ii=1,nR(iib,iparm)
 459 #ifdef MPI
 460               fimax_p(ii,iib,iparm)=v(1,ii,iib,iparm)+entfac(1)
 461               do t=2,scount(me)
 462                 if(v(t,ii,iib,iparm)+entfac(t).gt.fimax_p(ii,iib,iparm))
 463      &            fimax_p(ii,iib,iparm)=v(t,ii,iib,iparm)+entfac(t)
 464               enddo
 465 #else
 466               fimax(ii,iib,iparm)=v(1,ii,iib,iparm)+entfac(1)
 467               do t=2,ntot(islice)
 468                 if(v(t,ii,iib,iparm)+entfac(t).gt.fimax(ii,iib,iparm))
 469      &            fimax(ii,iib,iparm)=v(t,ii,iib,iparm)+entfac(t)
 470               enddo
 471 #endif
 472             enddo ! ii
 473           enddo ! iib
 474         enddo ! iparm
 475 #ifdef MPI
 476         call MPI_AllReduce(fimax_p(1,1,1),fimax(1,1,1),
 477      &   maxR*MaxT_h*nParmSet,MPI_DOUBLE_PRECISION,
 478      &   MPI_MAX,WHAM_COMM,IERROR)
 479 #endif
 480         do iparm=1,nParmSet
 481           do iib=1,nT_h(iparm)
 482             do ii=1,nR(iib,iparm)
 483 #ifdef MPI
 484               fi_p(ii,iib,iparm)=0.0d0
 485               do t=1,scount(me)
 486                 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
 487      &            +dexp(v(t,ii,iib,iparm)+entfac(t)-fimax(ii,iib,iparm))
 488 #ifdef DEBUG
 489               write (iout,'(4i5,3e15.5)') t,ii,iib,iparm,
 490      &         v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm)
 491 #endif
 492               enddo
 493 #else
 494               fi(ii,iib,iparm)=0.0d0
 495               do t=1,ntot(islice)
 496                 fi(ii,iib,iparm)=fi(ii,iib,iparm)
 497      &            +dexp(v(t,ii,iib,iparm)+entfac(t)-fimax(ii,iib,iparm))
 498               enddo
 499 #endif
 500             enddo ! ii
 501           enddo ! iib
 502         enddo ! iparm
 503
 504 #ifdef MPI
 505 #ifdef DEBUG
 506         write (iout,*) "fi before MPI_Reduce me",me,' master',master
 507         do iparm=1,nParmSet
 508           do ib=1,nT_h(nparmset)
 509             write (iout,*) "iparm",iparm," ib",ib
 510             write (iout,*) "beta=",beta_h(ib,iparm)
 511             write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
 512           enddo
 513         enddo
 514 #endif
 515 c        write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
 516 c     &   maxR*MaxT_h*nParmSet
 517 c        write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
 518 c     &   " WHAM_COMM",WHAM_COMM
 519         call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
 520      &   MPI_DOUBLE_PRECISION,
 521      &   MPI_SUM,Master,WHAM_COMM,IERROR)
 522 #ifdef DEBUG
 523         write (iout,*) "fi after MPI_Reduce nparmset",nparmset
 524         do iparm=1,nParmSet
 525           write (iout,*) "iparm",iparm
 526           do ib=1,nT_h(iparm)
 527             write (iout,*) "beta=",beta_h(ib,iparm)
 528             write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 529           enddo
 530         enddo
 531 #endif
 532         if (me1.eq.Master) then
 533 #endif
 534         avefi=0.0d0
 535         do iparm=1,nParmSet
 536           do ib=1,nT_h(iparm)
 537             do i=1,nR(ib,iparm)
 538               fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fimax(i,ib,iparm)
 539               avefi=avefi+fi(i,ib,iparm)
 540             enddo
 541           enddo
 542         enddo
 543         avefi=avefi/nFi
 544         do iparm=1,nParmSet
 545           write (iout,*) "Parameter set",iparm
 546           do ib =1,nT_h(iparm)
 547             write (iout,*) "beta=",beta_h(ib,iparm)
 548             do i=1,nR(ib,iparm)
 549               fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
 550             enddo
 551             write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 552             write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
 553           enddo
 554         enddo
 555
 556 ! Compute the norm of free-energy increments.
 557         finorm=0.0d0
 558         do iparm=1,nParmSet
 559           do ib=1,nT_h(iparm)
 560             do i=1,nR(ib,iparm)
 561               finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
 562               f(i,ib,iparm)=fi(i,ib,iparm)
 563             enddo
 564           enddo
 565         enddo
 566
 567         write (iout,*) 'Iteration',iter,' finorm',finorm
 568
 569 #ifdef MPI
 570         endif
 571         call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
 572      &   MPI_DOUBLE_PRECISION,Master,
 573      &   WHAM_COMM,IERROR)
 574         call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
 575      &   WHAM_COMM,IERROR)
 576 #endif
 577 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
 578         if (finorm.lt.fimin) then
 579           write (iout,*) 'Iteration converged'
 580           goto 20
 581         endif
 582
 583       enddo ! iter
 584
 585    20 continue
 586 ! Now, put together the histograms from all simulations, in order to get the
 587 ! unbiased total histogram.
 588 #ifdef MPI
 589       do t=0,tmax
 590         hfin_ent_p(t)=0.0d0
 591       enddo
 592 #else
 593       do t=0,tmax
 594         hfin_ent(t)=0.0d0
 595       enddo
 596 #endif
 597       write (iout,*) "--------------hist"
 598 #ifdef MPI
 599       do iparm=1,nParmSet
 600         do i=0,nGridT
 601           sumW_p(i,iparm)=0.0d0
 602           sumE_p(i,iparm)=0.0d0
 603           sumEbis_p(i,iparm)=0.0d0
 604           sumEsq_p(i,iparm)=0.0d0
 605           do j=1,nQ+2
 606             sumQ_p(j,i,iparm)=0.0d0
 607             sumQsq_p(j,i,iparm)=0.0d0
 608             sumEQ_p(j,i,iparm)=0.0d0
 609           enddo
 610         enddo
 611       enddo
 612       upindE_p=0
 613 #else
 614       do iparm=1,nParmSet
 615         do i=0,nGridT
 616           sumW(i,iparm)=0.0d0
 617           sumE(i,iparm)=0.0d0
 618           sumEbis(i,iparm)=0.0d0
 619           sumEsq(i,iparm)=0.0d0
 620           do j=1,nQ+2
 621             sumQ(j,i,iparm)=0.0d0
 622             sumQsq(j,i,iparm)=0.0d0
 623             sumEQ(j,i,iparm)=0.0d0
 624           enddo
 625         enddo
 626       enddo
 627       upindE=0
 628 #endif
 629 c 8/26/05 entropy distribution
 630 #ifdef MPI
 631       entmin_p=1.0d10
 632       entmax_p=-1.0d10
 633       do t=1,scount(me1)
 634 c        ent=-dlog(entfac(t))
 635         ent=entfac(t)
 636         if (ent.lt.entmin_p) entmin_p=ent
 637         if (ent.gt.entmax_p) entmax_p=ent
 638       enddo
 639       write (iout,*) "entmin",entmin_p," entmax",entmax_p
 640       call flush(iout)
 641       call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
 642      &  WHAM_COMM,IERROR)
 643       call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
 644      &  WHAM_COMM,IERROR)
 645       ientmax=entmax-entmin
 646       if (ientmax.gt.2000) ientmax=2000
 647       write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
 648       call flush(iout)
 649       do t=1,scount(me1)
 650 c        ient=-dlog(entfac(t))-entmin
 651         ient=entfac(t)-entmin
 652         if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
 653       enddo
 654       call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
 655      &  MPI_SUM,WHAM_COMM,IERROR)
 656       if (me1.eq.Master) then
 657         write (iout,*) "Entropy histogram"
 658         do i=0,ientmax
 659           write(iout,'(f15.4,i10)') entmin+i,histent(i)
 660         enddo
 661       endif
 662 #else
 663       entmin=1.0d10
 664       entmax=-1.0d10
 665       do t=1,ntot(islice)
 666         ent=entfac(t)
 667         if (ent.lt.entmin) entmin=ent
 668         if (ent.gt.entmax) entmax=ent
 669       enddo
 670       ientmax=-dlog(entmax)-entmin
 671       if (ientmax.gt.2000) ientmax=2000
 672       do t=1,ntot(islice)
 673         ient=entfac(t)-entmin
 674         if (ient.le.2000) histent(ient)=histent(ient)+1
 675       enddo
 676       write (iout,*) "Entropy histogram"
 677       do i=0,ientmax
 678         write(iout,'(2f15.4)') entmin+i,histent(i)
 679       enddo
 680 #endif
 681       do iparm=1,nParmSet
 682
 683         call restore_parm(iparm)
 684 c
 685 C Histograms
 686 c
 687 #ifdef MPI
 688         do ib=1,nT_h(iparm)
 689           do t=0,tmax
 690             hfin_p(t,ib)=0.0d0
 691           enddo
 692         enddo
 693         do i=1,maxindE
 694           histE_p(i)=0.0d0
 695         enddo
 696 #else
 697         do ib=1,nT_h(iparm)
 698           do t=0,tmax
 699             hfin(t,ib)=0.0d0
 700           enddo
 701         enddo
 702         do i=1,maxindE
 703           histE(i)=0.0d0
 704         enddo
 705 #endif
 706         do ib=1,nT_h(iparm)
 707           do i=0,MaxBinRms
 708             do j=0,MaxBinRgy
 709               hrmsrgy(j,i,ib)=0.0d0
 710 #ifdef MPI
 711               hrmsrgy_p(j,i,ib)=0.0d0
 712 #endif
 713             enddo
 714           enddo
 715         enddo
 716 #ifdef MPI
 717         do t=1,scount(me1)
 718 #else
 719         do t=1,ntot(islice)
 720 #endif
 721           indE = aint(potE(t,iparm)-aint(potEmin))
 722           if (indE.ge.0 .and. indE.le.maxinde) then
 723             if (indE.gt.upindE_p) upindE_p=indE
 724             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
 725           endif
 726 #ifdef MPI
 727           do ib=1,nT_h(iparm)
 728             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
 729             hfin_p(ind,ib)=hfin_p(ind,ib)+
 730      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
 731              if (rmsrgymap) then
 732                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
 733                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
 734                hrmsrgy_p(indrgy,indrms,ib)=
 735      &           hrmsrgy_p(indrgy,indrms,ib)+expfac
 736              endif
 737           enddo
 738 #else
 739           do ib=1,nT_h(iparm)
 740             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
 741             hfin(ind,ib)=hfin(ind,ib)+
 742      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
 743              if (rmsrgymap) then
 744                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
 745                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
 746                hrmsrgy(indrgy,indrms,ib)=
 747      &           hrmsrgy(indrgy,indrms,ib)+expfac
 748              endif
 749           enddo
 750 #endif
 751         enddo ! t
 752 c
 753 c Thermo and ensemble averages
 754 c
 755         do k=0,nGridT
 756           betaT=startGridT+k*delta_T
 757           call temp_scalfac(betaT,ft,ftprim,ftbis,*10)
 758 c            write (iout,*) "ftprim",ftprim
 759 c            write (iout,*) "ftbis",ftbis
 760           betaT=1.0d0/(1.987D-3*betaT)
 761 c 7/10/18 AL Determine the max Botzmann weights for each temerature
 762           call sum_ene(1,iparm,ft,etot)
 763           weimax(k,iparm)=-betaT*(etot-potEmin)+entfac(1)
 764 c          write (iout,*) "k",k," t",1," weight",weimax(k,iparm)
 765 #ifdef MPI
 766           do t=2,scount(me1)
 767 #else
 768           do t=2,ntot(islice)
 769 #endif
 770             call sum_ene(t,iparm,ft,etot)
 771             weight=-betaT*(etot-potEmin)+entfac(t)
 772 c            write (iout,*) "k",k," t",t," weight",weight
 773             if (weight.gt.weimax(k,iparm)) weimax(k,iparm)=weight
 774           enddo
 775 #ifdef MPI
 776         enddo
 777 #ifdef DEBUG
 778         write (iout,*) "weimax before REDUCE"
 779         write (iout,*) (weimax(k,iparm),k=0,ngridt)
 780 #endif
 781         do k=0,nGridT
 782           weimax_(k)=weimax(k,iparm)
 783         enddo
 784         call MPI_Allreduce(weimax_(0),weimax(0,iparm),nGridT+1,
 785      &   MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
 786 #ifdef DEBUG
 787         write (iout,*) "weimax"
 788         write (iout,*) (weimax(k,iparm),k=0,ngridt)
 789 #endif
 790         do k=0,nGridT
 791           temper=startGridT+k*delta_T
 792           betaT=1.0d0/(1.987D-3*temper)
 793           call temp_scalfac(temper,ft,ftprim,ftbis,*10)
 794           do t=1,scount(me1)
 795 #else
 796           do t=1,ntot(islice)
 797 #endif
 798             ind = ind_point(t)
 799 #ifdef MPI
 800             hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
 801 #else
 802             hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
 803 #endif
 804 c          write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
 805 c          call restore_parm(iparm)
 806             call sum_ene_deriv(t,iparm,ft,ftprim,ftbis,etot,eprim,ebis)
 807             weight=dexp(-betaT*(etot-potEmin)+entfac(t)-weimax(k,iparm))
 808 #ifdef DEBUG
 809             write (iout,*) "iparm",iparm," t",t," betaT",betaT,
 810      &        " etot",etot," entfac",entfac(t)," boltz",
 811      &        -betaT*(etot-potEmin)+entfac(t)," weimax",weimax(k,iparm),
 812      &        " weight",weight," ebis",ebis
 813 #endif
 814             etot=etot-temper*eprim
 815 #ifdef MPI
 816             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
 817             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
 818             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
 819             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
 820             do j=1,nQ+2
 821               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
 822               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
 823               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
 824      &         +etot*q(j,t)*weight
 825             enddo
 826 #else
 827             sumW(k,iparm)=sumW(k,iparm)+weight
 828             sumE(k,iparm)=sumE(k,iparm)+etot*weight
 829             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
 830             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
 831             do j=1,nQ+2
 832               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
 833               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
 834               sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
 835      &         +etot*q(j,t)*weight
 836             enddo
 837 #endif
 838           enddo ! t
 839         enddo ! k
 840         do ib=1,nT_h(iparm)
 841           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
 842      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
 843           if (rmsrgymap) then
 844           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
 845      &   (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
 846      &       WHAM_COMM,IERROR)
 847           endif
 848         enddo
 849         call MPI_Reduce(upindE_p,upindE,1,
 850      &    MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
 851         call MPI_Reduce(histE_p(0),histE(0),maxindE,
 852      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
 853
 854         if (me1.eq.master) then
 855
 856         if (histout) then
 857
 858         write (iout,'(6x,$)')
 859         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
 860      &   ib=1,nT_h(iparm))
 861         write (iout,*)
 862
 863         write (iout,'(/a)') 'Final histograms'
 864         if (histfile) then
 865           if (nslice.eq.1) then
 866             if (separate_parset) then
 867               write(licz3,"(bz,i3.3)") myparm
 868               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
 869             else
 870               histname=prefix(:ilen(prefix))//'.hist'
 871             endif
 872           else
 873             if (separate_parset) then
 874               write(licz3,"(bz,i3.3)") myparm
 875               histname=prefix(:ilen(prefix))//'_par'//licz3//
 876      &         '_slice_'//licz2//'.hist'
 877             else
 878               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
 879             endif
 880           endif
 881 #if defined(AIX) || defined(PGI)
 882           open (ihist,file=histname,position='append')
 883 #else
 884           open (ihist,file=histname,access='append')
 885 #endif
 886         endif
 887
 888         do t=0,tmax
 889           liczba=t
 890           sumH=0.0d0
 891           do ib=1,nT_h(iparm)
 892             sumH=sumH+hfin(t,ib)
 893           enddo
 894           if (sumH.gt.0.0d0) then
 895             do j=1,nQ
 896               jj = mod(liczba,nbin1)
 897               liczba=liczba/nbin1
 898               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
 899               if (histfile)
 900      &           write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
 901             enddo
 902             do ib=1,nT_h(iparm)
 903               write (iout,'(e20.10,$)') hfin(t,ib)
 904               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
 905             enddo
 906             write (iout,'(i5)') iparm
 907             if (histfile) write (ihist,'(i5)') iparm
 908           endif
 909         enddo
 910
 911         endif
 912
 913         if (entfile) then
 914           if (nslice.eq.1) then
 915             if (separate_parset) then
 916               write(licz3,"(bz,i3.3)") myparm
 917               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
 918             else
 919               histname=prefix(:ilen(prefix))//'.ent'
 920             endif
 921           else
 922             if (separate_parset) then
 923               write(licz3,"(bz,i3.3)") myparm
 924               histname=prefix(:ilen(prefix))//'par_'//licz3//
 925      &           '_slice_'//licz2//'.ent'
 926             else
 927               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
 928             endif
 929           endif
 930 #if defined(AIX) || defined(PGI)
 931           open (ihist,file=histname,position='append')
 932 #else
 933           open (ihist,file=histname,access='append')
 934 #endif
 935           write (ihist,'(a)') "# Microcanonical entropy"
 936           do i=0,upindE
 937             write (ihist,'(f8.0,$)') dint(potEmin)+i
 938             if (histE(i).gt.0.0e0) then
 939               write (ihist,'(f15.5,$)') dlog(histE(i))
 940             else
 941               write (ihist,'(f15.5,$)') 0.0d0
 942             endif
 943           enddo
 944           write (ihist,*)
 945           close(ihist)
 946         endif
 947         write (iout,*) "Microcanonical entropy"
 948         do i=0,upindE
 949           write (iout,'(f8.0,$)') dint(potEmin)+i
 950           if (histE(i).gt.0.0e0) then
 951             write (iout,'(f15.5,$)') dlog(histE(i))
 952           else
 953             write (iout,'(f15.5,$)') 0.0d0
 954           endif
 955           write (iout,*)
 956         enddo
 957         if (rmsrgymap) then
 958           if (nslice.eq.1) then
 959             if (separate_parset) then
 960               write(licz3,"(bz,i3.3)") myparm
 961               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
 962             else
 963               histname=prefix(:ilen(prefix))//'.rmsrgy'
 964             endif
 965           else
 966             if (separate_parset) then
 967               write(licz3,"(bz,i3.3)") myparm
 968               histname=prefix(:ilen(prefix))//'_par'//licz3//
 969      &         '_slice_'//licz2//'.rmsrgy'
 970             else
 971              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
 972             endif
 973           endif
 974 #if defined(AIX) || defined(PGI)
 975           open (ihist,file=histname,position='append')
 976 #else
 977           open (ihist,file=histname,access='append')
 978 #endif
 979           do i=0,nbin_rms
 980             do j=0,nbin_rgy
 981               write(ihist,'(2f8.2,$)')
 982      &          rgymin+deltrgy*j,rmsmin+deltrms*i
 983               do ib=1,nT_h(iparm)
 984                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
 985                   write(ihist,'(e14.5,$)')
 986      &              -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
 987      &              +potEmin
 988                 else
 989                   write(ihist,'(e14.5,$)') 1.0d6
 990                 endif
 991               enddo
 992               write (ihist,'(i2)') iparm
 993             enddo
 994           enddo
 995           close(ihist)
 996         endif
 997         endif
 998       enddo ! iparm
 999 #ifdef MPI
1000       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1001      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1002       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1003      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1004       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1005      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1006       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1007      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1008       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1009      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1010       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1011      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1012      &   WHAM_COMM,IERROR)
1013       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1014      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1015      &   WHAM_COMM,IERROR)
1016       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1017      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1018      &   WHAM_COMM,IERROR)
1019       if (me.eq.master) then
1020 #endif
1021       write (iout,'(/a)') 'Thermal characteristics of folding'
1022       if (nslice.eq.1) then
1023         nazwa=prefix
1024       else
1025         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1026       endif
1027       iln=ilen(nazwa)
1028       if (nparmset.eq.1 .and. .not.separate_parset) then
1029         nazwa=nazwa(:iln)//".thermal"
1030       else if (nparmset.eq.1 .and. separate_parset) then
1031         write(licz3,"(bz,i3.3)") myparm
1032         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1033       endif
1034       do iparm=1,nParmSet
1035       if (nparmset.gt.1) then
1036         write(licz3,"(bz,i3.3)") iparm
1037         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1038       endif
1039       open(34,file=nazwa)
1040       if (separate_parset) then
1041         write (iout,'(a,i3)') "Parameter set",myparm
1042       else
1043         write (iout,'(a,i3)') "Parameter set",iparm
1044       endif
1045       do i=0,NGridT
1046         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1047         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1048      &    sumW(i,iparm)
1049         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1050      &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1051         do j=1,nQ+2
1052           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1053           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1054      &     -sumQ(j,i,iparm)**2
1055           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1056      &     -sumQ(j,i,iparm)*sumE(i,iparm)
1057         enddo
1058         sumW(i,iparm)=(-dlog(sumW(i,iparm))-weimax(i,iparm))*(1.987D-3*
1059      &   (startGridT+i*delta_T))+potEmin
1060         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1061      &   sumW(i,iparm),sumE(i,iparm)
1062         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1063         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1064      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1065         write (iout,*)
1066         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1067      &   sumW(i,iparm),sumE(i,iparm)
1068         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1069         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1070      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1071         write (34,*)
1072         call flush(34)
1073       enddo
1074       close(34)
1075       enddo
1076       if (histout) then
1077       do t=0,tmax
1078         if (hfin_ent(t).gt.0.0d0) then
1079           liczba=t
1080           jj = mod(liczba,nbin1)
1081           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1082      &     hfin_ent(t)
1083           if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1084      &      dmin+(jj+0.5d0)*delta,
1085      &     hfin_ent(t)
1086         endif
1087       enddo
1088       if (histfile) close(ihist)
1089       endif
1090
1091 #ifdef ZSCORE
1092 ! Write data for zscore
1093       if (nslice.eq.1) then
1094         zscname=prefix(:ilen(prefix))//".zsc"
1095       else
1096         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1097       endif
1098 #if defined(AIX) || defined(PGI)
1099       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1100 #else
1101       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1102 #endif
1103       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1104       do iparm=1,nParmSet
1105         write (izsc,'("NT=",i1)') nT_h(iparm)
1106       do ib=1,nT_h(iparm)
1107         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1108      &    1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1109         jj = min0(nR(ib,iparm),7)
1110         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1111         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1112         write (izsc,'("&")')
1113         if (nR(ib,iparm).gt.7) then
1114           do ii=8,nR(ib,iparm),9
1115             jj = min0(nR(ib,iparm),ii+8)
1116             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1117             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1118             write (izsc,'("&")')
1119           enddo
1120         endif
1121         write (izsc,'("FI=",$)')
1122         jj=min0(nR(ib,iparm),7)
1123         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1124         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1125         write (izsc,'("&")')
1126         if (nR(ib,iparm).gt.7) then
1127           do ii=8,nR(ib,iparm),9
1128             jj = min0(nR(ib,iparm),ii+8)
1129             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1130             if (jj.eq.nR(ib,iparm)) then
1131               write (izsc,*)
1132             else
1133               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1134               write (izsc,'(t80,"&")')
1135             endif
1136           enddo
1137         endif
1138         do i=1,nR(ib,iparm)
1139           write (izsc,'("KH=",$)')
1140           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1141           write (izsc,'(" Q0=",$)')
1142           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1143           write (izsc,*)
1144         enddo
1145       enddo
1146       enddo
1147       close(izsc)
1148 #endif
1149 #ifdef MPI
1150       endif
1151 #endif
1152
1153       return
1154 C#undef DEBUG
1155    10 return1
1156       end
1157 c------------------------------------------------------------------------
1158       subroutine temp_scalfac(betaT,ft,ftprim,ftbis,*)
1159       implicit none
1160       include "DIMENSIONS"
1161       include "DIMENSIONS.FREE"
1162       include "COMMON.CONTROL"
1163       include "COMMON.FREE"
1164       include "COMMON.IOUNITS"
1165       double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
1166      &  kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
1167      &  logfac,tanhT,betaT,denom,eplus,eminus
1168       integer l
1169       if (rescale_mode.eq.1) then
1170         quot=betaT/T0
1171         quotl=1.0d0
1172         kfacl=1.0d0
1173         do l=1,5
1174           quotl1=quotl
1175           quotl=quotl*quot
1176           kfacl=kfacl*kfac
1177           denom=kfacl-1.0d0+quotl
1178           fT(l)=kfacl/denom
1179           ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
1180           ftbis(l)=l*kfacl*quotl1*
1181      &     (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
1182         enddo
1183 #if defined(FUNCTH)
1184         ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1185      &            320.0d0
1186         ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1187         ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1188      &          /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1189 #elif defined(FUNCT)
1190         fT(6)=betaT/T0
1191         ftprim(6)=1.0d0/T0
1192         ftbis(6)=0.0d0
1193 #else
1194         fT(6)=1.0d0
1195         ftprim(6)=0.0d0
1196         ftbis(6)=0.0d0
1197 #endif
1198       else if (rescale_mode.eq.2) then
1199         quot=betaT/T0
1200         quotl=1.0d0
1201         do l=1,5
1202           quotl1=quotl
1203           quotl=quotl*quot
1204           eplus=dexp(quotl)
1205           eminus=dexp(-quotl)
1206           logfac=1.0d0/dlog(eplus+eminus)
1207           tanhT=(eplus-eminus)/(eplus+eminus)
1208           fT(l)=1.12692801104297249644d0*logfac
1209           ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1210           ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
1211      &    2*l*quotl1/T0*logfac*
1212      &    (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
1213      &    +ftprim(l)*tanhT)
1214         enddo
1215 #if defined(FUNCTH)
1216         ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1217      &           320.0d0
1218         ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1219         ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1220      &         /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1221 #elif defined(FUNCT)
1222         fT(6)=betaT/T0
1223         ftprim(6)=1.0d0/T0
1224         ftbis(6)=0.0d0
1225 #else
1226         fT(6)=1.0d0
1227         ftprim(6)=0.0d0
1228         ftbis(6)=0.0d0
1229 #endif
1230       else if (rescale_mode.eq.0) then
1231         do l=1,5
1232           fT(l)=1.0d0
1233           ftprim(l)=0.0d0
1234         enddo
1235       else
1236         write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
1237      &    rescale_mode
1238         call flush(iout)
1239         return1
1240       endif
1241       return
1242       end
1243 c--------------------------------------------------------------------
1244       subroutine sum_ene(t,iparm,ft,etot)
1245       implicit none
1246       include 'DIMENSIONS'
1247       include 'DIMENSIONS.ZSCOPT'
1248       include 'DIMENSIONS.FREE'
1249       include 'COMMON.CONTROL'
1250       include 'COMMON.FFIELD'
1251       include "COMMON.SBRIDGE"
1252       include "COMMON.ENERGIES"
1253       include "COMMON.IOUNITS"
1254       integer t,iparm
1255       double precision fT(6)
1256       double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
1257      &  escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
1258      &  eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
1259      &  eliptran,esaxs
1260       evdw=enetb(21,t,iparm)
1261       evdw_t=enetb(1,t,iparm)
1262 #ifdef SCP14
1263       evdw2_14=enetb(17,t,iparm)
1264       evdw2=enetb(2,t,iparm)+evdw2_14
1265 #else
1266       evdw2=enetb(2,t,iparm)
1267       evdw2_14=0.0d0
1268 #endif
1269 #ifdef SPLITELE
1270       ees=enetb(3,t,iparm)
1271       evdw1=enetb(16,t,iparm)
1272 #else
1273       ees=enetb(3,t,iparm)
1274       evdw1=0.0d0
1275 #endif
1276       ecorr=enetb(4,t,iparm)
1277       ecorr5=enetb(5,t,iparm)
1278       ecorr6=enetb(6,t,iparm)
1279       eel_loc=enetb(7,t,iparm)
1280       eello_turn3=enetb(8,t,iparm)
1281       eello_turn4=enetb(9,t,iparm)
1282       eturn6=enetb(10,t,iparm)
1283       ebe=enetb(11,t,iparm)
1284       escloc=enetb(12,t,iparm)
1285       etors=enetb(13,t,iparm)
1286       etors_d=enetb(14,t,iparm)
1287       ehpb=enetb(15,t,iparm)
1288       estr=enetb(18,t,iparm)
1289       esccor=enetb(19,t,iparm)
1290       edihcnstr=enetb(20,t,iparm)
1291       eliptran=enetb(22,t,iparm)
1292       esaxs=enetb(26,t,iparm)
1293 #ifdef SPLITELE
1294       if (shield_mode.gt.0) then
1295         etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1296      &    +ft(1)*welec*ees
1297      &    +ft(1)*wvdwpp*evdw1
1298      &    +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1299 c     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1300      &    +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1301      &    +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1302      &    +ft(2)*wturn3*eello_turn3
1303      &    +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1304      &    +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1305      &    +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1306       else
1307         etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1308      &    +wvdwpp*evdw1
1309      &    +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1310 c     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1311      &    +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1312      &    +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1313      &    +ft(2)*wturn3*eello_turn3
1314      &    +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1315      &    +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1316      &    +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1317       endif
1318 #else
1319       if (shield_mode.gt.0) then
1320         etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1321      &    +ft(1)*welec*(ees+evdw1)
1322      &    +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1323 c     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1324      &    +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1325      &    +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1326      &    +ft(2)*wturn3*eello_turn3
1327      &    +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1328      &    +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1329      &    +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1330       else
1331         etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1332      &    +ft(1)*welec*(ees+evdw1)
1333      &    +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1334 c     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1335      &    +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1336      &    +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1337      &    +ft(2)*wturn3*eello_turn3
1338      &    +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1339      &    +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1340      &    +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1341       endif
1342 #endif
1343       return
1344       end
1345 c--------------------------------------------------------------------
1346       subroutine sum_ene_deriv(t,iparm,ft,ftprim,ftbis,etot,eprim,ebis)
1347       implicit none
1348       include 'DIMENSIONS'
1349       include 'DIMENSIONS.ZSCOPT'
1350       include 'DIMENSIONS.FREE'
1351       include 'COMMON.CONTROL'
1352       include 'COMMON.FFIELD'
1353       include "COMMON.SBRIDGE"
1354       include 'COMMON.ENERGIES'
1355       include "COMMON.IOUNITS"
1356       integer t,iparm
1357       double precision fT(6),fTprim(6),fTbis(6),
1358      &  eprim,ebis,temper
1359       double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
1360      &  escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
1361      &  eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
1362      &  eliptran,esaxs
1363       evdw=enetb(21,t,iparm)
1364       evdw_t=enetb(1,t,iparm)
1365 #ifdef SCP14
1366       evdw2_14=enetb(17,t,iparm)
1367       evdw2=enetb(2,t,iparm)+evdw2_14
1368 #else
1369       evdw2=enetb(2,t,iparm)
1370       evdw2_14=0.0d0
1371 #endif
1372 #ifdef SPLITELE
1373       ees=enetb(3,t,iparm)
1374       evdw1=enetb(16,t,iparm)
1375 #else
1376       ees=enetb(3,t,iparm)
1377       evdw1=0.0d0
1378 #endif
1379       ecorr=enetb(4,t,iparm)
1380       ecorr5=enetb(5,t,iparm)
1381       ecorr6=enetb(6,t,iparm)
1382       eel_loc=enetb(7,t,iparm)
1383       eello_turn3=enetb(8,t,iparm)
1384       eello_turn4=enetb(9,t,iparm)
1385       eturn6=enetb(10,t,iparm)
1386       ebe=enetb(11,t,iparm)
1387       escloc=enetb(12,t,iparm)
1388       etors=enetb(13,t,iparm)
1389       etors_d=enetb(14,t,iparm)
1390       ehpb=enetb(15,t,iparm)
1391       estr=enetb(18,t,iparm)
1392       esccor=enetb(19,t,iparm)
1393       edihcnstr=enetb(20,t,iparm)
1394       eliptran=enetb(22,t,iparm)
1395       esaxs=enetb(26,t,iparm)
1396 #ifdef SPLITELE
1397       if (shield_mode.gt.0) then
1398         etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1399      &    +ft(1)*welec*ees
1400      &    +ft(1)*wvdwpp*evdw1
1401      &    +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1402 c     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1403      &    +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1404      &    +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1405      &    +ft(2)*wturn3*eello_turn3
1406      &    +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1407      &    +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1408      &    +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1409         eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
1410 C     &            +ftprim(6)*evdw_t
1411      &    +ftprim(1)*wscp*evdw2
1412      &    +ftprim(1)*welec*ees
1413      &    +ftprim(1)*wvdwpp*evdw1
1414      &    +ftprim(1)*wtor*etors+
1415      &    ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1416      &    ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1417      &    ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1418      &    ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1419      &    ftprim(1)*wsccor*esccor
1420         ebis=ftbis(1)*wsc*(evdw+ft(6)*evdw_t)
1421      &    +ftbis(1)*wscp*evdw2+
1422      &    ftbis(1)*welec*ees
1423      &    +ftbis(1)*wvdwpp*evdw
1424      &    +ftbis(1)*wtor*etors+
1425      &    ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1426      &    ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1427      &    ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1428      &    ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1429      &    ftbis(1)*wsccor*esccor
1430       else
1431         etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1432      &    +wvdwpp*evdw1
1433      &    +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1434 c     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1435      &    +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1436      &    +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1437      &    +ft(2)*wturn3*eello_turn3
1438      &    +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1439      &    +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1440      &    +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1441         eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
1442      &    +ftprim(1)*wtor*etors+
1443      &    ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1444      &    ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1445      &    ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1446      &    ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1447      &    ftprim(1)*wsccor*esccor
1448         ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
1449      &    ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1450      &    ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1451      &    ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1452      &    ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1453      &    ftbis(1)*wsccor*esccor
1454       endif
1455 #else
1456       if (shield_mode.gt.0) then
1457         etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1458      &    +ft(1)*welec*(ees+evdw1)
1459      &    +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1460 c     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1461      &    +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1462      &    +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1463      &    +ft(2)*wturn3*eello_turn3
1464      &    +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1465      &    +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1466      &    +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1467         eprim=ftprim(1)*(evdw+ft(6)*evdw_t)
1468      &    +ftprim(1)*welec*(ees+evdw1)
1469      &    +ftprim(1)*wtor*etors+
1470      &     ftprim(1)*wscp*evdw2+
1471      &     ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1472      &     ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1473      &     ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1474      &     ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1475      &     ftprim(1)*wsccor*esccor
1476        ebis= ftbis(1)*(evdw+ft(6)*evdw_t)
1477      &    +ftbis(1)*wscp*evdw2
1478      &    +ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1479      &     ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1480      &     ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1481      &     ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1482      &     ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1483      &     ftprim(1)*wsccor*esccor
1484       else
1485         etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1486      &    +ft(1)*welec*(ees+evdw1)
1487      &    +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1488 c     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1489      &    +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1490      &    +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1491      &    +ft(2)*wturn3*eello_turn3
1492      &    +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1493      &    +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1494      &    +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1495         eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
1496      &    +ftprim(1)*wtor*etors+
1497      &     ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1498      &     ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1499      &     ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1500      &     ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1501      &     ftprim(1)*wsccor*esccor
1502         ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1503      &     ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1504      &     ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1505      &     ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1506      &     ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1507      &     ftprim(1)*wsccor*esccor
1508       endif
1509 #endif
1510       return
1511       end