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