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