Ask Question

Name:
Title:
Your Question:

Answer Question

Name:
Your Answer:
User Submitted Source Code!


Description:
  bds
Language: FORTRAN
Code:
     Module Bspglocons ! Bead Spring simulation, Global constant (parameters)
  Save  ! always use to ensure contents remein unchanged between calls

  !-- Physical
  Integer, Parameter :: Ndim = 3 ! Dimension of simulation
  Integer, Parameter ::  MaxXdist = 101

  !-- Computational dimensions
  Integer, Parameter :: k4b = Selected_int_kind(9)
  Integer, Parameter :: SNGL = Selected_real_kind(4)
  Integer, Parameter :: DOBL = Selected_real_kind(8)

  ! Choose the precision of computation here
  !Integer, Parameter :: FPREC = DOBL
  Integer, Parameter :: FPREC = SNGL

  !-- Math Constants
  Real (FPREC) , Parameter :: PI = 3.14159265358979323846
  Real (FPREC) , Parameter :: TINI = 1e-25
  Real (FPREC) , Parameter :: MYEPS = 1e-6
  Integer, Parameter :: MAXCHEB = 500


  Integer, Parameter :: Hook = 1, FENE = 2, ILC = 3, WLC = 4

  !-- Averages
  Integer, Parameter :: NAvgs = 63
  ! for Eqbm and Prod runs of TrajType
  Integer, Parameter :: EQ = 1, PR = 2

  !-- Time correlation
  Integer, Parameter :: Ncorrel = 4
  Integer, Parameter :: ALLOCAT=0, INIT=1, SAMPLE=2, PRNT=3  &
       , GETDAT=4, SAVDAT=5

  !-- for ran_1()
  Integer, Parameter :: RESET=1, SHOW=2

End Module Bspglocons


Module Bspglovars ! Bead Spring simulation, Global variables
  Use Bspglocons
  Save  ! always use to ensure contents remein unchanged between calls




  Real (FPREC)  Imploop_tol
End Module Bspglovars


Module BlockTcorr
  Use Bspglocons
  Save  ! always use to ensure contents remain unchanged between calls

  Integer, Parameter :: MAXBLKCOR = 2 ! Maximum number of blocks
  Integer NTOPB 

  Real (FPREC) , Parameter :: TDIFFMAX = 3e6

  Integer   Nto(MAXBLKCOR)     &     ! Number of elements filled
       , Iblmc                       ! counter to track MaxBlock
  

  Integer, Allocatable :: Npp(:,:)   ! number of TC sampled
  Real (FPREC), Allocatable ::   PPcorr(:,:,:) ! correl funcs 

End Module BlockTcorr



Module Bspintfacs
  Use Bspglocons

  !---------------------------------------------------------------------c     
  !     Driver subroutines                                              c
  !---------------------------------------------------------------------c

  Interface
     Subroutine Initial_position(Stype,N,L0s,R,seed)
       Use bspglocons
       Integer, Intent (in) :: Stype
       Integer(k4b), Intent (in) :: N
       Real (FPREC) , Intent (in) :: L0s
       Integer(k4b), Intent (inout) :: seed
       Real (FPREC) , Intent (out) :: R(:,:)
       !Real, intent (out) :: R(Ndim,N)
     End Subroutine Initial_position
  End Interface

  Interface
     Subroutine GaussRand(N,GX,seed)
       Use bspglocons
       Integer(k4b), Intent (in) :: N
       Integer(k4b), Intent (inout) :: seed
       Real (FPREC) , Intent (out), Dimension(N) :: GX
     End Subroutine GaussRand
  End Interface


  Interface
     Subroutine Thermalise_Chain(NBeads, R_Bead, spring_type, &
          Ntcur, Ntsteps, Delts, &
          Hstar,chbjerrium,ldebye,zstar,rcstar,dstar,cutoffL, L0s, FHI, &
          TrajType, nseed, Ntcsamp, Ntsamp, propsums, Nsdone,ifindex)
       Use bspglocons
       Implicit None
       Integer, Intent(in) :: NBeads
       Real (FPREC) , Intent (inout), Dimension(:,:) :: R_Bead 

       Integer, Intent (in) :: spring_type ! Spring force law type
       ! Hookean = 1, FENE = 2, ILC = 3, WLC = 4

       Integer(k4b) , Intent (in) :: Ntcur ! Current time steps completed
       Integer , Intent (in) :: Ntsteps ! time steps to be advanced to
       Real (FPREC) , Intent (in) :: Delts  ! Integration time interval specs


       Real (FPREC) ,  Intent (in) :: Hstar         ! HI parameter (non-dim)
! EV parameters (non-dim)
       Real (FPREC) ,  Intent (in) :: zstar, rcstar, dstar,cutoffL
! Electrostatic parameters (non-dim)
       Real (FPREC) ,  Intent (in) :: chbjerrium, ldebye
       Real (FPREC) ,  Intent (in) :: L0s           ! finite ext., param sqrt(b)
       Logical , Intent (in) :: FHI  ! fluctuating/preaveraged HI 

       Integer, Intent (in) :: Trajtype ! EQuilibriation, PRoduction

       Integer (k4b), Intent (inout) :: nseed 

       Integer, Intent (in), Optional :: Ntcsamp, Ntsamp,ifindex
       Real (FPREC), Intent (inout), Optional :: propsums(:) 
       Integer, Intent (inout), Optional :: Nsdone 

     End Subroutine Thermalise_Chain


  End Interface


  Interface
     Subroutine polish_poly_root(c,xin,atol)
       Use Bspglocons
       Implicit None
       ! use newton raphson to polish the root of a polynomial
       Integer, Parameter :: n=4
       Real (FPREC) , Intent (in) :: c(n)
       Real (FPREC) , Intent (inout) :: xin
       Real (FPREC) , Intent (in) :: atol
     End Subroutine polish_poly_root
  End Interface

  !---------------------------------------------------------------------c     
  !     Utility subroutines                                             c
  !---------------------------------------------------------------------c

  Interface
     Subroutine ran_1(action, n, R, idum)
       Use bspglocons
       Integer, Intent(in) :: action
       Integer(k4b), Intent(in) :: n
       Real (FPREC) , Intent(inout) :: R(n)
       Integer(k4b), Intent(inout) ::idum
     End Subroutine ran_1
  End Interface

  Interface
     Subroutine maxminev_fi(n, A, maxev, minev)
       Use Bspglocons
       Integer n
       Real (FPREC)  A(:,:,:,:), maxev, minev
     End Subroutine maxminev_fi
  End Interface

  Interface
     Subroutine chbyshv (L, d_a, d_b, a)
       Use bspglocons
       Integer L
       Real (FPREC)  d_a, d_b, a(0:MAXCHEB)
     End Subroutine chbyshv
  End Interface

  Interface
     Subroutine cublu
       Use bspglocons
     End Subroutine cublu
  End Interface
  !---------------------------------------------------------------------c
  !     Properties estimators                                           c
  !---------------------------------------------------------------------c     

  Interface
     Subroutine eqchain_props(N, Rbead, F, propsums, Rg)
       Use bspglocons
       Implicit None
       Integer, Intent (in) :: N
       Real (FPREC) , Intent (in), Dimension(:,:)  :: Rbead, F
       Real (FPREC) , Intent (inout),Dimension(:) ::propsums, Rg(Ndim)
     End Subroutine eqchain_props
  End Interface

  Interface
     Subroutine numint(f,dx,nmin,nmax,nord,sumf)
       Use Bspglocons
       Implicit None

       Real (FPREC) , Intent(in), Dimension(:) :: f
       Real (FPREC) , Intent(in) :: dx
       Real (FPREC) , Intent(out) :: sumf
       Integer, Intent (in) :: nmin, nmax,nord
     End Subroutine numint
  End Interface

  Interface
     Subroutine meanerr(vec,mean,err)
       Use Bspglocons
       Real (FPREC) , Intent(in) :: vec(:)
       Real (FPREC) , Intent(out) :: mean
       Real (FPREC) , Intent(out) :: err
     End Subroutine meanerr
  End Interface

  Interface
     Function normeqpr(Stype, r, b)
       Use Bspglocons
       Implicit None
       Integer, Intent (in) :: Stype
       Real (FPREC) , Intent (in) :: b, r
       Real (FPREC)   normeqpr
     End Function normeqpr
  End Interface

  Interface
     Function  lam1_th(hs, N)
       Use Bspglocons
       Real (FPREC)  lam1_th
       Real (FPREC) ,  Intent (in) :: hs
       Integer,  Intent (in) :: N
     End Function lam1_th
  End Interface



  Interface
     Subroutine Sample_Tcorr(action, sfile, pval)
       Use Bspglocons
       Use BlockTcorr
       Implicit None

       Integer, Intent(in) :: action
       Real (FPREC) , Intent(in), Optional :: pval(Ndim,Ncorrel) 
       Character (*), Intent (in), Optional :: sfile

     End Subroutine Sample_Tcorr

  End Interface

End Module Bspintfacs

Program Eqbm2Ensemble

!!!====================   Declarations ===========================!!!

  !! *** Load Modules containing global variables and interfaces *** !!

  Use bspglocons
  Use bspglovars
  !Use BlockSums
  Use BlockTcorr
  Use bspintfacs
  Implicit None

  !!*** Variables  ***!!

!!$  !*** Additional debugging information
!!$  !Logical, parameter :: DEBUG = .TRUE. 
!!$  Logical, Parameter :: DEBUG = .False. 
!!$  Integer, Parameter :: debunit=16
!!$  Real (8) wtime, wtime0, wtimeinit


  !*** parallel processing                  
 ! Include 'mpif.h'
  Integer mpierr, Nprocs, MyId


  !*** Model parameters and variables
  Integer  :: Stype, NBeads
  Real (FPREC) :: hstar, sqrtb
! EV parameters (non-dim)
       Real (FPREC) :: zstar, rcstar, dstar,cutoffL
! Electrostatic parameters (non-dim)
       Real (FPREC) :: chbjerrium, ldebye
  Real (FPREC), Allocatable :: PosVecR(:,:), pv_main(:,:), pv_cv(:,:)

  !*** Time-Integration  parameters

  !--
  !  Set the multiple of Rouse relaxation times used for the 
  !  EQuilibriation trajectory, one with equilibrium time step
  !  and another with production run time step
  !--

  Integer, parameter :: EqEqTime = 8, EqPrTime = 1
  Real (FPREC), parameter :: deltseq = 1

  Integer :: Ntraj, Ntrtot, Ntsteps
  Integer (k4b) :: Nitot, Ntdone 
  Real (FPREC) :: deltspr


  !*** Sampling parameters
  Integer  :: Maxlam, TrajType, Ntsamp, Ntcsamp, Nsdone
  Real (FPREC) :: Epslam, Dctime

  !*** File i/o
  Integer, parameter :: SavMins = 15 ! time interval between saves
  Character (len=20) infile, outfile, gavgfile, corfile, gppfile &
       , resfile, contfile, xcfile
  Integer, Parameter ::   stdout=5, inunit=10, resunit=11, corunit = 12 &
       , gavunit=13, gppunit=14, cunit=17, xcunit=19, avunit=20
  Character (len=10), Parameter :: ErString="Error"


  !*** Trajectories related
  Logical Filexists, ToContinue, FHI, Flocked
  Integer (k4b) nseed, nseed_cv
  Real (FPREC) ::  Delts, Conf_t, dumX(1)

  !*** Misc
  Character(10) :: Prop_names(NAvgs), Correl_names(Ncorrel)
  Integer i, j, nprops,continueprog

  Real (FPREC) tlongest,  t1rouse, t1zimm

  Real (FPREC), Allocatable :: avgs(:)

  Integer ni,  iblock, nblock, gntrajdone, gntrthis, ntrajproc
  Real (FPREC), Allocatable, Dimension(:,:,:) :: xcorravgs,  gxcavgs
  Real (FPREC), Allocatable, Dimension(:,:) :: Xvar, Yvar


  Integer ppcount, ib,ifindex
  Real (FPREC) ctime  
  Real (FPREC), Allocatable ::  ppbuf(:,:,:), Gppavg(:,:,:,:) 

  Real time_begin, time_sa, time_se, time_end


!!!=================== Execution begin ================================!!!

!  !!*** MPI Intialisation ***!!
!
!  Call MPI_Init(mpierr)                    
!
!  !-- Obtain the total number of processors
!  Call MPI_Comm_Size (MPI_Comm_World, Nprocs, mpierr) 
!
!  !-- Obtain the rank or ID of the current processor
!  Call MPI_Comm_Rank (MPI_Comm_World, MyId, mpierr)                

   Nprocs = 1
   MyId = 0


  !!*** Input DATA ***!!
  infile="par-eq.in"
  Open (unit=inunit,file=infile,status="old")
  Read (inunit,*)
  Read (inunit,*) SType, NBeads, hstar, zstar, rcstar, sqrtb
  Read (inunit,*)
  Read (inunit,*) Ntraj, Ntrtot, Nitot, deltspr, Imploop_tol
  !^^ Note Imploop_tol is a global variable
  Read (inunit,*)
  Read (inunit,*) Epslam, Maxlam, dstar, chbjerrium, ldebye,cutoffL,continueprog
  Close (unit=inunit)



  !!*** Memory Allocation  ***!!

  !-- longest relaxation times obtained from Rouse and Zimm Models
  t1rouse =   0.5/Sin(PI/2/Nbeads)**2
  If (hstar >= 0.52) Then
  t1zimm = lam1_th(0.5,NBeads) ! Thurstons approximation
  else
  t1zimm = lam1_th(hstar,NBeads) ! Thurstons approximation
  End If

  !-- Number of time steps per sample (for time correl)
  Ntcsamp = Max(t1zimm * epslam/deltspr, 1._FPREC) 

  !--Number of time origins per block (number of points per block)
  NTOPB = Maxlam * t1zimm / (Ntcsamp * deltspr)

  If (NTOPB <= 1) Then
     print *, 'Error: NTOPB=', NTOPB, ' is <=1. Check the input'
     !call MPI_Abort(mpierr)
     stop
  End If
  



  Allocate(PosVecR(Ndim,NBeads), pv_main(Ndim,NBeads), pv_cv(Ndim,NBeads) )
  Allocate ( &
       avgs(NAvgs)  &
       )

  Allocate ( &
       xcorravgs(NAvgs,Nitot,3) &
       , Xvar(Navgs,Nitot) &
       , Yvar(Navgs,Nitot) &
       )

  

  call Sample_Tcorr(ALLOCAT)



  !*** Variables specific to root processor only

  If (MyId.Eq.0) Then

     Call CPU_Time(time_begin)

     Allocate ( &
          gxcavgs(NAvgs,Nitot,3)  &
          !-- receive buffer
          , ppbuf(Ncorrel,MAXBLKCOR,0:NTOPB) &  
          !-- global avg, for main and cv
          , Gppavg(Ncorrel,MAXBLKCOR,0:NTOPB,2) & 
          )

     ! the standard error of mean obtained from t-distribution
     ! for sample mean and sample standard deviation
     Conf_t = 2 ! 95% confidence for degrees of freedom > 20

     Prop_names(1) = "R^2"        ! square of end-to-end vector          
     Prop_names(2) = "Rg^2"       ! sq. Radius of gyration
     Prop_names(4) = "Rinv"       ! Sum <1/rmn> / N^2
     Prop_names(5) = "TrDiff"  ! Center of mass diffusivity
      Prop_names(3) = "Lp"  ! persistence length
     Prop_names(6) = "X11"          ! Avg Stretch 
     Prop_names(7) = "X12"  ! Avg Stretchx2
     Prop_names(8) = "X13"  !Avg Stretch  x3
     Prop_names(9) = "X22"  ! structure factor
     Prop_names(10) = "X23"  ! structure factor
     Prop_names(11) = "X33"  ! structure factor
     do i=12,NAvgs
      Prop_names(i) = "Sq"
     enddo

     Correl_names(1) = "Re"
     Correl_names(2) = "A=D.F/4"
     Correl_names(3) = "Sxy"

  End If


  !!*** Initialisation ***!!

  !*** Generate a pseudo random seed
 if(continueprog.eq.1)then  
call GetaSeed(nseed)
else
open(unit=97,file='posveccont',status='old')
read(97,*)nseed
close(97)
endif

  !--debug suntemp nseed = 201271
  print *, 'Initial seed of ', MyId, ' = ', nseed 



  !*** Other initialisations

  !-- the current number of time steps completed (not used in this 
  !-- algorithm)
  Ntdone = 0

  !-- Filenames
  gavgfile = "gavgseq.dat"
  gppfile  = "gppcorr.dat"

  resfile = "result.dat" ! only for the Master processor 
  contfile = "continue"
  xcfile = 'xcorr.dat' ! cross correlation function
  corfile = "acorr.dat"


  !*** Continuation run (if possible)

  !-- obtain the # completed trajectories from the disk
  Inquire (file=gavgfile, exist=Filexists)
  If (Filexists) Then
     Open (unit=gavunit,file=gavgfile,status='old')
     Read (gavunit,*) gntrajdone
     Close(gavunit)
  Else
     gntrajdone = 0
  End If

  !-- Ntraj is the number of trajectories for this run
  !-- Ntrtot is the total trajs desired
  Ntraj = Min(Ntraj,Ntrtot-gntrajdone)

  !-- trajectories per processor
  nblock = Max(Ntraj/Nprocs,1) ! Note integer division



  !-- Initialisation configuration of the main variate's PosVec
     if(continueprog.eq.1)then
  Call Initial_position(SType,Nbeads,sqrtb,pv_main,nseed)
     else
     open(unit=95,file='posvecread',status='old')
     read(95,*)pv_main
     close(95)
     endif


  !-- 
  !  Number of time steps: for one sampling production run,
  !  and for sampling time correlation function
  !--
  Ntsteps = t1zimm/deltspr ! approx decorrelation time steps
  Ntsamp = Ntsteps ! sample once after decorrelation

  !-- sampling interval for time-correl
  Dctime = deltspr*Ntcsamp
  



  !-- averages for the two variates and their cross correlation function
  xcorravgs = 0
  !-- number of trajectories completed in the current run (in each processor)
  ntrajproc = 0


  !-- number of items to be transferred in mpi_reduce
  ppcount = Ncorrel*MAXBLKCOR*(NTOPB+1)

  !-- Global time correlation fn data
  If (MyID ==0) Gppavg = 0

  !-- Start time to save data at regular intervals
  Call CPU_Time(time_sa)



     !--
     PosVecR = pv_main   

     !-- save for future use for control variate below
     pv_cv = PosVecR  

     !-- 
     !  for a variation across trajectories,'cos ran_1 is reset below.
     !  It is possible that nseed remains unaltered across one trajectory,
     !  if sufficient number of calls are not made.
     !--
!     nseed = nseed+1  
     nseed_cv = nseed
     Call ran_1(RESET,1,dumX,nseed)


     !*** Equilibriate chain configuration 

     !--
     !  An exact synchronisation of the control variate with the main is
     !  not valid for cases other than Hookean in theta solvent.  In all 
     !  other cases, it is essential to equilibriate the control variate
     !  to its equilibrium distribution, starting from the synchronised value.
     !  The sampling is done only after this.  In order to do this both 
     !  the main and the CV are thermalised for a predetermined duration
     !  of time (typically about 10 T_Zimm).  Since we have seen
     !  earlier that the variates remain correlated for more than 300 
     !  T_Zimm, we can be sure that after 10 Tzimm, both will remain 
     !  correlated, and obey the respective equilibrium distribution 
     !  functions.  Note the equilibriation is required only for the 
     !  control variate, however, in order to maintain the same RN 
     !  sequence, the main also needs to go through the same sequence.
     !--



!     If (zstar .Ne. 0 .or. SType /= Hook) Then
!        TrajType = EQ

        !--debug print *, 'eqbting in main'
        
!        Ntsteps = EqEqTime*t1rouse/deltseq

        !-- Free draining at a larger time step

!        Call Thermalise_Chain(NBeads, PosVecR, SType,  &
!             Ntdone, Ntsteps, deltspr, &
!             0._FPREC,chbjerrium,ldebye,zstar,rcstar,dstar,cutoffL,&
!      sqrtb, FHI, TrajType, nseed, &
!             Ntcsamp, Ntsamp, avgs, Nsdone,ifindex)

        !-- Free draining at the production time step
!        Ntsteps = EqPrTime*t1rouse/deltspr

!        Call Thermalise_Chain(NBeads, PosVecR, SType,  &
!             Ntdone, Ntsteps, deltspr, &
!             0._FPREC,chbjerrium,ldebye,zstar,rcstar,dstar,cutoffL,&
!      sqrtb, FHI, TrajType, nseed, &
!             Ntcsamp, Ntsamp, avgs, Nsdone,ifindex)
        
!     End If

     !--debug Call ran_1(SHOW,1,dumX,nseed)
     !--debug print *, MyID, 'RNG status (main):', dumX
     
     !***  Production run

     FHI=.True.

     !-- Production trajectory for sampling
     TrajType = PR
     
     !-- Number of time steps to thermalise
     Ntsteps = t1zimm/deltspr ! approx decorrelation time
     Ntsamp = Ntsteps ! one sample after decorrelation

     !-- Initialise static averages and time corr variables
     Xvar = 0
     avgs = 0
     Call Sample_Tcorr(INIT) ! PPcorr is initialised to 0

     Do ni=1,Nitot
     ifindex=ni
!     Xvar(:,ni)=avgs 
        !--
        !  avgs is a cumulative sum in the following call; to get the 
        !  instantaeneous value it is added here and will be 
        !  subtracted out below. Note that Ntsamp = Ntsteps implies
        !  that only one sample is taken for one call.(at the begining)
        !--

        !-- Actual Brownian dynamics done only here!
   !  Ntsteps = 10.0*t1zimm/deltspr ! approx decorrelation time
        Call Thermalise_Chain(NBeads, PosVecR, SType,  &
             Ntdone, Ntsteps, deltspr, &
             hstar,chbjerrium,ldebye,zstar,rcstar,dstar,cutoffL,&
      sqrtb, FHI, TrajType, nseed, &
             Ntcsamp, Ntsamp, avgs, Nsdone,ifindex)

     Open (unit=96,file='posvecstore',status="unknown")
     write(96,*)PosVecR
     close(96)
!        Xvar(:,ni) = avgs - Xvar(:,ni)
        Xvar(:,ni) = avgs 
     End Do

     !-- Time correlations from sample_tcorr
     Do i=1,Ncorrel
        Where (Npp > 0) &
             PPcorr(i,:,:) = PPcorr(i,:,:) / (Npp / Ncorrel)
           ! since Npp is incremented by Ncorrel for each sample
     End Do

     !Call MPI_Barrier (MPI_Comm_World, mpierr)
     !Call GMPI_Redsum(PPcorr,ppbuf,ppcount,mpierr)
     ppbuf = PPCorr


     !-- 
     !  transfer to the global average of main as ppbuf will be
     !  overwritten 
     !--
     If (MyID == 0) Gppavg(:,:,:,1) = Gppavg(:,:,:,1) + ppbuf


     !-- 
     !  save for the next trajectory. Note: PosVecR will b 
     !  overwritten in the control variate below
     !--
     pv_main = PosVecR 



     !-- store for the mean and cross-correlation 
     xcorravgs(:,:,1) =  xcorravgs(:,:,1) +  Xvar
 !     xcorravgs(:,:,2) =  xcorravgs(:,:,2) +  Yvar
 !    xcorravgs(:,:,3) =  xcorravgs(:,:,3) +  Xvar*Yvar

     !-- number of trajectories completed in this processor
     ntrajproc = ntrajproc + 1

     !-- 
     !  While it is expected that for equilibrium runs, atleast,
     !  each of the processors will take the same amount of time 
     !  irrespective of the starting configuration,  we however impose a
     !  barrier here to ensure synchronised computation.  This is because
     !  it is possible (in principle) that one of the processors 
     !  completes all the 'nblock' trajectories within the time 
     !  saving interval and exits the loop.  This would mean that the
     !  consolidation of data across processors done in 'xcorrsav()' 
     !  will be from different calls, one following and one outside the loop.
     !  The "fast" processors will end up calling xcorrsav only once!
     !--
     !Call MPI_Barrier (MPI_Comm_World, mpierr)
     !--debug print *, 'END trajectory : ', iblock, Myid

     !--End time for writing to disk at regular intervals
     Call CPU_Time(time_se)

     !-- Save data at regular intervals
     If ( (time_se-time_sa)/60  >= SavMins) Then
     !If ( (time_se-time_sa)  >= 5) Then
        !--
        !  the internal procedure xcorrsav, adds the current data
        !  to that in the disk (if present). Note: in the process
        !  the variables xcorravgs, ntrajproc is overwritten and 
        !  is reset to zero in the end.  Since xcorrsav involves
        !  consolidation of data across processors, no specific
        !  processor condition is used 
        !--

        !-- 
        !  no of trajdone cud vary across procs for the same time interval,
        !  a simple multiple of Nprocs will not hold good in general
        !  Note that since ntrajproc is transferred, it will be used
        !  as a buffer in the following calls; it nolonger retains
        !  the meaning of its name!
        !--
        !Call MPI_REDUCE (ntrajproc, gntrthis, 1,      &
        !     MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, mpierr)
        gntrthis = ntrajproc

        Call xcorrsav

        Call PPcorrsav


        time_sa  = time_se
     End If !-- to save or not to save


  !*** Final consolidation
  !Call MPI_REDUCE (ntrajproc, gntrthis, 1,      &
  !     MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, mpierr)
  gntrthis = ntrajproc

  Call xcorrsav
  Call PPcorrsav


  If (MyID == 0) Then
     Call CPU_Time(time_end)
     Print *, 'CPU_time usage:',  (time_end-time_begin)/3600.0, 'hrs'
  End If

  
  !--debug print *, Myid, iblock
  !Call MPI_Finalize(mpierr)      


!!! ========================== Format Statements ========================== !!!

1 Format ('#',16(A10,2X))
  !2  Format ('#', 2(I10,2X), 12(G10.3,2X), (I10,2X)  )
2 Format ('#', 2(G2.0,10X), 12(G11.4,1X), (G10.3,2X), G2.0,11X  )

3 Format ('#', <2*NAvgs+4>(A11,3X))
  ! short cut way to get (nearly) left justification
31 Format ('#', <2*NAvgs+4>(G2.0,12X))

4 Format (<2*NAvgs+4>(G13.6,1X))


72 Format (5(F10.4,2x))
80 Format (3(F10.4,2x))


200 Format ('#',6(G12.0,2X),<nprops>(G12.0,2X, 'Error       ',2X))
202 Format ('#',<6+nprops*2>(I2,12X)) !approx LJustification
210 Format (I3,  11X, <6+nprops*2>(G12.5,2X))

  !250  Format (3(2x,f13.6))
250 Format (3(2x,g13.6))


300 Format(G11.4,2x,A6,G11.4,A25)
301 Format(G11.4,2x,A6,G11.4,A25,2x, I15)



Contains 

!!! =======================================================================
!!!                    Internal Procedures Begin
!!! =======================================================================

Subroutine GetaSeed(aseed)
  Integer (k4b) aseed
  Character (len=12) clk(3)
  Integer clok(8)
  
  !-- clock(8) gives the milliseconds
  Call Date_and_time(clk(1), clk(2), clk(3), clok) 
  aseed = clok(8)*100000+clok(7)*1000+clok(6)*10+clok(5)

  !-- an unique seed incase clok is the same
  aseed = (aseed + 201271)*(MyId +1)  

  !-- shuffling, lifted from ran_1
  aseed = Ieor(aseed, Ishft(aseed,13))
  aseed = Ieor(aseed, Ishft(aseed,-17))
  aseed = Ieor(aseed, Ishft(aseed,5))

End Subroutine GetaSeed



!  Subroutine GMPI_Redsum(var,buf,N,mpierr)
!    !--
!    !  This procedure is a generic wrapper for MPI_Reduce(),
!    !  for various floating point precisions
!    !--
!    Use Bspglocons
!    Implicit None
!
!    Integer N,mpierr 
!    Real (FPREC):: var(N), buf(N)
!
!    Include 'mpif.h'
!
!    !--debug print *, 'count requested = ', N, '. size(var) = ', size(var)
!    if (N /= size(var)) Then
!       print *, 'Mismatch of requested size and buffer size', N, size(var)
!       return
!    end if
!    
!
!    Select Case(FPREC)
!    Case(SNGL)
!       Call MPI_Reduce (var, buf, N,      &
!            MPI_Real, MPI_Sum, 0, MPI_Comm_World, mpierr)
!    Case(DOBL)
!       Call MPI_Reduce (var, buf, N,      &
!            MPI_Double_Precision, MPI_Sum, 0, MPI_Comm_World, mpierr)
!    Case DEFAULT
!       Write (*,*) 'Wrong precision in GMPI_redsum()'
!       Call MPI_Abort(MPI_Comm_World,mpierr)
!    End Select
! End Subroutine GMPI_Redsum



  Subroutine xcorrsav
    Integer count

    !-- -----------
    !  This internal procedure, adds the current static data to that 
    !  in the disk (if present), thru the following steps
    !  * consolidate current results from processors 
    !       xcorravgs -> gxcavgs
    !  * read data from disk into xcorravgs
    !  * add to gxcavgs
    !  * write gxcavgs to disk
    !  * reset xcorravgs to zero
    !
    !  Note on usage:  
    !    Since a consolidation step is involved, this procedure must
    !    be called by all processors, and not exclusive to a particular MyID.
    !-- -----------

    !!*** consolidate results from the other processors *** !!

    count = Size(xcorravgs)

    !Call GMPI_Redsum(xcorravgs,gxcavgs,count,mpierr)
    gxcavgs = xcorravgs


    If (MyID == 0) Then

       !-- obtain global values from the disk
       Inquire (file=gavgfile, exist=Filexists)
       If (Filexists) Then
          !-- check if it is locked
          !Flocktest: Do 
          !   Call islocked(gavgfile,Flocked)
           !  If(.Not. Flocked) Then
           !     Call lockfile(gavunit,gavgfile) 
                !^^ to be unlocked only after new data is written

                Open (unit=gavunit,file=gavgfile,status='old')

                !-- 
                !  read in ntrajproc,xcorravgs as the current run's 
                !  values have already been transferred to global values
                !--
                Read (gavunit,*) ntrajproc
                Read (gavunit,*) xcorravgs
                Close (gavunit)
           !     Exit Flocktest
           !  Else  ! when File is locked
                ! Sleep for a while and try again until it is unlocked,  
                ! requires -Vaxlib in ifort
           !     Call sleep(1)
           !  End If ! file lock test

         ! End Do Flocktest
       Else ! if Globalavgs File_does_not_exist
          !-- lock the file b4 opening a new one for writing
         ! Call lockfile(gavunit,gavgfile) 
          !^^ tobe unlocked after new data is written

          ntrajproc = 0
          xcorravgs = 0
       End If

       !-- Note that the 2nd term in RHS is obtained from the disk
       gxcavgs = gxcavgs + xcorravgs
       gntrajdone = gntrthis + ntrajproc
       !-- gntrthis was consolidated in the main program, b4 this call


       !--
       !  Save the updated global averages data on to the disk
       !  Note that the gavunit file is still locked till this
       !  is completed
       !--

       Open(unit=gavunit,file=gavgfile,status='replace')

       Write (gavunit,*) gntrajdone
       Write (gavunit,*) gxcavgs

       Write (gavunit,*) 'End-of-File'
       Close(gavunit)

       !-- Release global avgs file for use to other programs (jobs)
     !  Call unlockfile(gavunit,gavgfile)


       !-- Compute averages and cross correlations
       gxcavgs = gxcavgs/gntrajdone

       Open (unit = xcunit, file = xcfile, status="unknown")
       Write (xcunit, 590) '# Cross correlation functions separated one t1zimm.'
       Write (xcunit, 592) '# Trajectories: ', gntrajdone
       Write (xcunit, 600) "# Time     ", (Prop_names(i), i = 1,NAvgs-1) 

       Do ni = 1,Nitot
          Write (xcunit, 601) ni, &
               ( (gxcavgs(i,ni,3) -  &
               gxcavgs(i,ni,1) * gxcavgs(i,ni,2)), &
               i = 1,Navgs-1)
       End Do
       Close(xcunit)


       Open (unit = avunit, file = "avgs.dat", status="unknown")
Write(avunit, 640) '#Ensemble averages of functions by one t1zimm, for the FHI and PREAV.' 
       Write (avunit,642) '# Trajectories:', gntrajdone

       Write (avunit, 650 ) "# Time       ", (Prop_names(i), Prop_names(i) &
            , i = 1,NAvgs-1) 

       Do ni = 1,Nitot
          Write (avunit, 651) ni, &
               (  gxcavgs(i,ni,1) , gxcavgs(i,ni,2) , &
               i = 1,Navgs-1)
       End Do
       Close(avunit)

    End If ! Root processor


    !-- 
    !  Since the data has been written to disk the following
    !  variables need to be reset before continuing with the simulation
    !--
    xcorravgs = 0
    ntrajproc = 0



590 Format (A)
592 Format (A,G12.5)
600 Format (A, 2x, <NAvgs>(A13,2x))
601 Format (G13.6, 2x, <NAvgs>(G13.6, 2x) )

640 Format (A)
642 Format (A,G12.5)
650 Format (A, 2x, <2*(NAvgs)>(A13,2x))
651 Format (G13.6, 2x, <2*(NAvgs)>(G13.6, 2x) )

  End Subroutine xcorrsav


  Subroutine PPcorrsav

    !-- -----------
    !  This internal procedure, adds the current Time correl data to that 
    !  in the disk (if present), thru the following steps
    !  * read data from disk into ppbuf
    !  * add to gppavg
    !  * write gppavg to disk
    !  * write time correl to disk
    !  * reset gppavg to zero
    !
    !-- -----------

    !-- 
    !  Note that in the case of PPcorr, the results have been consolidated
    !  from each processor after each trajectory
    !-- 



    If (MyID ==0) Then

       !-- obtain global values from the disk
       Inquire (file=gppfile, exist=Filexists)
       If (Filexists) Then
          !-- check if it is locked
          !Flocktest: Do 
            ! Call islocked(gppfile,Flocked)
            ! If(.Not. Flocked) Then
            !    Call lockfile(gppunit,gppfile) 
                !^^ to be unlocked only after new data is written

                Open (unit=gppunit,file=gppfile,status='old')

                !-- 
                !  read in ntrajproc, ppbuf as the current run's 
                !  values have already been transferred to global avgs
                !--

                Read (gppunit,*) ntrajproc
                gntrajdone = gntrthis + ntrajproc
                !-- gntrthis was consolidated in the main program, b4 this call

                !-- 
                !  Given that the total length and interval of
                !  sampling remains same for main and cv, and
                !  across several runs, the maximum sampled
                !  points in a block Nto(ib) (see sample_tcorr),
                !  remains same.  So it does not matter when it is
                !  read into the same variable.  However, not reading
                !  this will not print out the time correlation function
                !  in corfile for the case: when no trajectories
                !  are sampled in this run and there is only data in
                !  the disk to be processed.
                !--
                Read (gppunit,*) Iblmc, Nto

                Read (gppunit,*) ppbuf
                Gppavg(:,:,:,1) = Gppavg(:,:,:,1) + ppbuf

                !-- for the CV
                Read (gppunit,*) ppbuf
                Gppavg(:,:,:,2) = Gppavg(:,:,:,2) + ppbuf

                Close (gppunit)
          !      Exit Flocktest
          !   Else  ! when File is locked
                ! Sleep for a while and try again until it is unlocked,  
                ! requires -Vaxlib in ifort
           !     Call sleep(1)
           !  End If ! file lock test

         ! End Do Flocktest
       Else ! if Globalavgs File_does_not_exist
          !-- lock the file b4 opening a new one for writing
         ! Call lockfile(gppunit,gppfile) 
          !^^ tobe unlocked after new data is written

          !-- These statements are redundant, but just there to
          !-- show the similarities with xcorrsav()
          !-- ntrajproc = 0
          !-- ppbuf = 0 
       End If ! gppfile exists

       !--
       !  Save the new global averages data on to the disk
       !  Note that the gppunit file is still locked till this
       !  is completed
       !--

       Open(unit=gppunit,file=gppfile,status='replace')

       Write (gppunit,*) gntrajdone
       Write (gppunit,*) Iblmc, Nto
       Write (gppunit,*) Gppavg(:,:,:,1)
       Write (gppunit,*) Gppavg(:,:,:,2)

       Write (gppunit,*) 'End-of-File'
       Close(gppunit)

       !-- Release global avgs file for use to other programs (jobs)
       !Call unlockfile(gppunit,gavgfile)


       !**** Calculate time correlation functions

       Gppavg = Gppavg/gntrajdone

       !--
       !  The stress correlation function comes from six components
       !  of the stress tensor, and their dot product is stored
       !  in Ncorrel=3 and Ncorrel=4.  Average them out here,
       !  and store in Ncorrel=3
       !--
       Gppavg(3,:,:,:) = (Gppavg(3,:,:,:) + Gppavg(3,:,:,:))/6

       Open (unit = corunit, file = corfile, status="unknown")
       Write (corunit,*) '# Time correlation functions for the main and CV'
       Write (corunit,*) '# Trajectories completed:', gntrajdone
       Write (corunit,500) "Time       ", "Time/t1zimm"  &
            , (correl_names(i), correl_names(i), i = 1,Ncorrel-1)

       ib=1
       Do j = 0, 1
          ctime = j*Dctime*(NTOPB**(ib-1))
          Write (corunit, 501 ) ctime, ctime/t1zimm &
               , (Gppavg(i,ib,j,1), Gppavg(i,ib,j,2), i = 1,Ncorrel-1)
       End Do

       Do ib = 1, Min(MAXBLKCOR, Iblmc)
          Do j = 2, Min(Nto(ib), NTOPB)
             ctime = j*Dctime*(NTOPB**(ib-1))
             Write (corunit, 501 ) ctime, ctime/t1zimm &
                  , (Gppavg(i,ib,j,1), Gppavg(i,ib,j,2), i = 1,Ncorrel-1)
          End Do
       End Do

       Close(corunit)

       !!*** Write continuation information to disk ***!!

       contfile="continue"
       Open(unit=cunit,file=contfile,status='replace')

       If (Ntrtot - gntrajdone > 0) Then
          Write (cunit,*) (Ntrtot - gntrajdone) &
               , '  more trajectories to be completed of ', Ntrtot
          Close(cunit)
       Else 
          Close(cunit,status='delete')
       End If

       !-- Note that gppavg is defined only in the root proc
       Gppavg = 0

    End If ! Root processor


    !-- 
    !  Since the data has been written to disk the following
    !  variables need to be reset before continuing with the simulation
    !  Gppavg, howver needs to b reset only in the root processor.
    !  Setting it to zero in other procs (without allocation)
    !  leads to the process hanging
    !--
    ntrajproc = 0



500 Format ('#', 2(A13, 2x), <2*(Ncorrel-1)>(A13, 2x))
501 Format (2(G13.6, 2x), <2*(Ncorrel-1)>(G13.6, 2x) )

  End Subroutine PPcorrsav


!!! =================================================================== !!!
!!!                    Internal Procedures End                          !!!
!!! =================================================================== !!!


End Program Eqbm2Ensemble



! some conventions followed for comments

!!!========== Main divisions (Parts) ============!!!

!!***   Section  ***!!
!*** Subsection

!-- Comment for the line(s) below
!^^ Comment for the line(s) above

!--
!   Block comment paragraph 
!   spanning
!   several 
!   lines 
!--

!!! Time-stamp: <utils.f90 10:51, 29 Apr 2005 by P Sunthar>

!_______________________________________
!  Contains general purpose utilities
!_______________________________________

!!! $Log: utils.f90,v $
!!! Revision 1.2  2005/05/24 01:07:10  pxs565
!!! variable NTOPB
!!!
!!! Revision 1.1  2005/03/16 00:54:06  pxs565
!!! Initial revision
!!!



Subroutine Initial_position(SType,N,L0s,R,myseed)
  Use bspglocons
  Use bspglovars
  Use bspintfacs
  !Use iflport ! for portable intel intrinsic functions, here ran()
  Implicit None

  Integer, Intent (in) :: Stype
  Integer(k4b), Intent (in) :: N
  Integer(k4b), Intent (inout) :: myseed
  Real (FPREC) , Intent (in) :: L0s
  Real (FPREC) , Intent (out), Dimension(:,:) :: R

  Integer, Save :: varseed = 201235

  Integer mu, done
  Real (FPREC) , Parameter :: upbound = 2
  Real (FPREC)  modr,b2b(Ndim),rv,eqpr, rvchk, b

  Call GaussRand(Ndim*N,R,myseed)



  R(:,1) = 0

  ! intrinsic ran() modifies the seed for every call
  ! so use a separate varible seed, so that the ran_1 sequence is
  ! not altered by this variation.
  !varseed = seed 

  b = L0s*L0s


!!$  do rv = 0,0.99,.01
!!$     eqpr = normeqpr(Stype,rv,b);
!!$     write (*,*) rv,eqpr
!!$  end do



  Do mu = 2,N
     b2b = R(:,mu) ! * sqrt(sigma=1), gaussian dist betwn two adj beads


     If (SType .Ne. HOOK) Then

        done = 0;
        Do While (done .Ne. 1)

           !rv = ran(varseed)
        call Random_Number(rv) ! standards conforming IBM
           eqpr = normeqpr(Stype,rv,b);


           ! use another uniform random number to decide acceptance
           ! ie, accept only eqpr fraction at this r
           !rvchk = ran(varseed)
        call Random_Number(rvchk) ! standards conforming IBM
           !write (*,*) rv, rvchk , eqpr
           If (rvchk .Le. eqpr/upbound) done = 1
        End Do

        modr = Sqrt(Sum(b2b*b2b))
        b2b = b2b/modr * (rv * L0s)
     End If

     R(:,mu) = R(:,mu-1) + b2b 

     ! in the case of a finitely extensible spring, a gaussian distribution
     ! can lead to incorrect forces some exceptional b2b vectors, we limit
     ! this by altering the distance when preseving the random direction.
     ! and an additional random factor between 0 and 1
     ! A good guess of the factor cud b obtained more rigorously.
  End Do
End Subroutine Initial_position

Subroutine GaussRand(N,GX,seed)
!!! Returns a normalised gaussian random variable with mean zero 
!!! and variance 1: f(r) = exp(-r^2/2) / sqrt(2*Pi)
!!! To obtain a dist with a given sigma use:  sqrt(sigma) * Gx

  Use bspglocons
  Use bspintfacs
  Implicit None

  Integer(k4b), Intent (in) :: N
  Integer(k4b), Intent (inout) :: seed
  Real (FPREC) , Intent (out), Dimension(N) :: GX

  Real (FPREC)  rnd(2), x,y, rsq, fac
  Integer i,n2
  Real (FPREC)  :: rtemp(N+1) 


  Do i = 1,N,2
     rsq = 0.0
     Do While (rsq .Ge. 1.0  .Or.  rsq .Eq. 0.0)
        Call ran_1(0,2,rnd,seed)
        x = 2*rnd(1) -1
        y = 2*rnd(2) -1
        rsq = x*x + y*y
     End Do
     fac = Sqrt(-2*Log(rsq)/rsq);
     rtemp(i) = x*fac
     rtemp(i+1) = y*fac
  End Do

  GX = rtemp(1:N)

End Subroutine GaussRand




Subroutine ran_1(action, n, X, idum)
  Use bspglocons
  Implicit None
  Integer, Intent(in) :: action
  Integer(k4b), Intent(inout) ::idum
  Integer(k4b), Intent(in) :: n
  Real (FPREC) , Intent(inout) :: X(n)
  !---------------------------------------------------------------------c
  !     Random number generator based on procedure given in             c
  !     Numerical Recipes in Fortran 90, Chapter B7, pp. 1141-1143      c
  !     This generator is supposed to have a period of about 3.1E18.    c
  !                                                                     c
  !     The calling sequence is the same as that of RANL,               c 
  !     the BLAS random number generator. The arguments are:            c
  !     n - gives the dimension of vector R                             c
  !     X - 1D array of dimension n; on exit                            c
  !         contains n unifromly distributed random numbers in (0,1)    c
  !     idum |- seeds which are updated on exit                         c
  !---------------------------------------------------------------------c

  Integer(k4b),Parameter :: IA=16807,IM=2147483647,IQ=127773,IR=2836
  Real (FPREC) , Save :: am
  Integer(k4b), Save :: ix=1, iy = -1, k, callcntr=0
  Integer(k4b) count


  Select Case (action)
  Case (RESET)
     ix = 1
     iy = -1
     k = 0
     callcntr = 0

  Case (SHOW)
     !-- pass the value of ix, which can be used to display
     X(1) = ix
     !-- Write (*,*) 'ran_1() status: ', callcntr,ix,iy,idum


  Case DEFAULT
     callcntr=callcntr+1
     If ((idum.Le.0) .Or. (iy.Lt.0)) Then
        am = Nearest(1.0,-1.0)/IM
        iy = Ior(Ieor(888889999,Abs(idum)),1)
        ix = Ieor(777755555,Abs(idum))
        idum = Abs(idum) + 1
     End If

     Do count = 1,n
        ix = Ieor(ix,Ishft(ix,13))
        ix = Ieor(ix, Ishft(ix,-17))
        ix = Ieor(ix,Ishft(ix, 5))
        k = iy/IQ
        iy = IA*(iy-k*IQ)-IR*k
        If (iy.Lt.0) iy = iy + IM
        X(count) = am*Ior(Iand(IM,Ieor(ix,iy)),1)
     End Do

  End Select

End Subroutine ran_1


Subroutine chbyshv (L, d_a, d_b, a)
  Use bspglocons
  Implicit None
  Integer L
  Real (FPREC)  d_a, d_b, a(0:MAXCHEB)
  !---------------------------------------------------------------------c
  !     This routine calculates the Chebyshev coefficients for the      c
  !     Chebyshev polynomial approximation of a square root function.   c
  !     The routine computes L+1 coefficients, which are returned in    c
  !     in the one-dimensional array a, for the approximation of the    c
  !     square root of any variable that lies in the interval           c 
  !     (d_a, d_b).                                                     c
  !---------------------------------------------------------------------c 


  Integer j, k
  Real (FPREC)  xks(0:L)


  !     Calculate the shift factors
  !     d_a = 2/(lmax-lmin)
  !      -d_b/d_a = (lmax+lmin)/2 

  !     Calculate the collocation points
  Do k = 0,L
     xks(k) = Cos(PI*(k+0.5)/(L+1))/d_a -d_b/d_a
  End Do

  !     Calculate the Chebyshev coefficients
  a = 0.0
  Do j = 0,L
     Do k = 0,L
        a(j) = a(j)+(xks(k)**0.5)*Cos(j*(k+0.5)*PI/(L+1))
     End Do
     a(j) = 2.0/(L+1)*a(j)
  End Do
  a(0) = a(0)/2

End Subroutine chbyshv



Subroutine maxminev_fi(n, A, maxev, minev)
  Use bspglocons
  Implicit None
  Integer n
  Real (FPREC)  A(:,:,:,:), maxev, minev
  !---------------------------------------------------------------------c
  !     This routine approximates the maximum and minimum eigen values  c
  !     of a symmetric matrix nxn matrix A using Fixman's suggestion.   c
  !---------------------------------------------------------------------c
  Integer lda, incx, incy, i
  Real (FPREC)  F(n), G(n), alpha, beta
  Real (SNGL) sdot
  Real (DOBL) ddot !! Blas

  F = 1.0
  G = 0.0

  lda = n
  alpha = 1.0
  beta = 0.0
  incx = 1
  incy = 1

  Select Case (FPREC)
  Case(SNGL)
     Call ssymv('U',n,alpha,A,lda,F,incx,beta,G,incy)
     maxev = sdot(n, F, 1, G, 1)
  Case(DOBL)
     Call dsymv('U',n,alpha,A,lda,F,incx,beta,G,incy)
     maxev = ddot(n, F, 1, G, 1)
  End Select


  maxev = 2.0*maxev/n

  !Forall (i = 1:n) F(i) = (-1.0)**i
  Do i = 1,n,2
     F(i)   = -1.
     F(i+1) =  1.
  End Do


  Select Case (FPREC)
  Case(SNGL)
     Call ssymv('U',n,alpha,A,lda,F,incx,beta,G,incy)
     minev = sdot(n, F, 1, G, 1)
  Case(DOBL)
     Call dsymv('U',n,alpha,A,lda,F,incx,beta,G,incy)
     minev = ddot(n, F, 1, G, 1)
  End Select



  minev = minev/2.0/n
  !        write (*,*) maxev, minev

End Subroutine maxminev_fi




Subroutine polish_poly_root(c,xin,atol)
  Use Bspglocons
  Implicit None
  ! use newton raphson to polish the root of a polynomial
  Integer, Parameter :: n=4 ! presently only for cubic
  Real (FPREC) , Intent (in) :: c(n)
  Real (FPREC) , Intent (inout) :: xin
  Real (FPREC) , Intent (in) :: atol

  Real (FPREC)  ::  p, p1,x0,x

  Integer i,iter
  Integer, Parameter :: IMAX = 30

  x = xin

  ! algo from NR: techniques for polishing, roots of polynomials
  Do iter=1,IMAX
     x0 = x
     p =  c(n)*x + c(n-1)
     p1 = c(n)
     Do i=n-2,1,-1
        p1 = p + p1*x
        p  = c(i) + p*x
     End Do

     !if (abs(p1) < atol) then ! f' = 0 is not solvable in newt-raph
     If (Abs(p1) .Eq. 0.0) Then ! f' = 0 is not solvable in newt-raph
        !! for WLC, in this case f'' also = 0, therefore
        p1 = 6 * c(4) ! f''' 
        p1 = 6*p/p1
        ! note sign(a,b) returns sgn(b) * mod(a)
        x = x - Sign(1._FPREC,p1) * Abs(p1)**(1./3.)
        !! we omit considering cases for ILC and FENE, as they
        !! donot seem to have this singularity
     Else
        x = x - p/p1
     End If

     If (Abs(p) .Le. atol) Exit
  End Do
  If (iter>IMAX) Then
     !write (5,*) 'poly-root: loop exceeded'
  End If


  xin = x
End Subroutine polish_poly_root


Subroutine numint(f,dx,nmin,nmax,nord,sumf)
  Use Bspglocons
  Implicit None

  Real (FPREC) , Intent(in), Dimension(:) :: f
  Real (FPREC) , Intent(in) :: dx
  Real (FPREC) , Intent(out) :: sumf
  Integer, Intent (in) :: nmin, nmax,nord



  Integer, Parameter :: stdout=5, stderr=6
  Integer nbound,nint,j

  nbound = Ubound(f,1)
  nint = nmax-nmin

  !write (*,*) 'nbound = ', nbound


  If (nint > nbound .Or. nmin < 1 .Or. nmax > nbound ) Then
     Write(stderr,*) 'Array out of bounds: ', nmin,nmax,nbound
     Return
  End If

  sumf = 0.

  Select Case (nord) 
  Case (1) ! trapezoidal rule
     sumf = 0.5 * (f(nmax) + f(nmin))

  Case(2)
     sumf = 5./12 * (f(nmin) + f(nmax)) + &
          13./12 *(f(nmin+1) + f(nmax-1)) 
  Case (3)
     sumf = 3./8 * (f(nmin) + f(nmax)) + &
          7./6 *(f(nmin+1) + f(nmax-1)) + &
          23./24 *(f(nmin+2) + f(nmax-2))
  End Select

  Do j = nmin+nord, nmax-nord
     sumf = sumf + f(j)
  End Do
  sumf =  dx*sumf

End Subroutine numint

Subroutine meanerr(vec,mean,err)
  Use Bspglocons
  Real (FPREC), Intent(in) :: vec(:)
  Real (FPREC), Intent(out) :: mean
  Real (FPREC), Intent(out) :: err ! standard error of mean

  Integer n

  n = Ubound(vec,1)
  mean = Sum(vec)/n

  ! though the following is correct, 
  ! err = sqrt((sum(vec*vec) - mean*mean*n)/(n-1))
  ! it does not always numerically evaluate to a quantity > 0 
  ! for very small errors, such as in the case of diffusivity.  
  ! so use the sure shot formula.
  err = 2*Sqrt(Sum((vec-mean)*(vec-mean))/(n-1)/n)
  ! the standard error of mean is approximately in the 95% confidence
  ! interval, giving the factor 2 (from the t-distribution)


End Subroutine meanerr


Function normeqpr(Stype, r, b)
  Use bspglocons
  Implicit None
  Integer, Intent (in) :: Stype
  Real (FPREC) , Intent (in) :: b, r
  Real (FPREC)   normeqpr

  Real (FPREC)  expphi
  External expphi

  Integer, Save :: norm_computed = 0
  Real (FPREC)    , Save :: norm_b, bsav

  Integer , Parameter :: Ln = 21 ! 21 always see below for X,W data
  Real (DOBL) X(Ln), W(Ln), scratch(Ln)
  Real (DOBL)  endpts(2)

  Integer n
  Real (FPREC)  xfact,rg, fb , M

  If ( norm_computed .Eq. 0  .And. bsav .Ne. b) Then

     !! generate the gauss abscissa and weights for legendre
     !call gaussq(1, Ln, 0, 0, 0, endpts, scratch, X, W) ;

     !! generated from c code
     X(01) = -0.993752170620389;   W(01) = 0.016017228257774
     X(02) = -0.967226838566306;   W(02) = 0.03695378977085309
     X(03) = -0.920099334150401;   W(03) = 0.05713442542685689
     X(04) = -0.853363364583318;   W(04) = 0.0761001136283793
     X(05) = -0.768439963475678;   W(05) = 0.09344442345603385
     X(06) = -0.667138804197413;   W(06) = 0.1087972991671478
     X(07) = -0.55161883588722;   W(07) = 0.1218314160537286
     X(08) = -0.424342120207439;   W(08) = 0.1322689386333376
     X(09) = -0.288021316802401;   W(09) = 0.1398873947910733
     X(10) = -0.145561854160895;   W(10) = 0.14452440398997
     X(11) = -2.4782829604619e-16;   W(11) = 0.1460811336496907
     X(12) = 0.145561854160895;   W(12) = 0.1445244039899697
     X(13) = 0.288021316802401;   W(13) = 0.1398873947910732
     X(14) = 0.424342120207439;   W(14) = 0.1322689386333371
     X(15) = 0.55161883588722;   W(15) = 0.1218314160537288
     X(16) = 0.667138804197413;   W(16) = 0.1087972991671492
     X(17) = 0.768439963475678;   W(17) = 0.09344442345603389
     X(18) = 0.853363364583317;   W(18) = 0.07610011362837897
     X(19) = 0.920099334150401;   W(19) = 0.05713442542685797
     X(20) = 0.967226838566306;   W(20) = 0.03695378977085323
     X(21) = 0.993752170620389;   W(21) = 0.01601722825777395




     M = 18  ! some large number, obtained by trial and error

     ! limit the upper bound of M
     If (M > b/2) M = b/2

     xfact = Sqrt(2*M/b)

     norm_b = 0
     fb = 0
     Do n=1,Ln
        rg = (X(n) + 1)/2 * xfact
        norm_b = norm_b +  W(n) * expphi(Stype,rg,b) * rg * rg
        fb = fb +  W(n) * expphi(Stype,rg,b) * rg**4
     End Do

     !Write (*,*) '# b = ', b
     !write (*,*) '# computed normb = ', norm_b
     !Write (*,*) '# f(b) = ', b*fb/norm_b

     bsav = b
     norm_computed = 1
  End If

  normeqpr = r*r*expphi(Stype,r,b)/norm_b

End Function normeqpr

Function expphi(Stype, r, b)
!!! return the distribution function (without normalisation)
  Use bspglocons
  Implicit None
  Integer, Intent (in) :: Stype
  Real (FPREC) , Intent (in) :: r, b
  Real (FPREC)   expphi

  Select Case (Stype)
  Case (HOOK)
     Write (*,*) 'Hookean case called in expb: Check the code'
  Case (FENE)
     expphi = (1-r*r)**(b/2)

  case (ILC)
     expphi =(exp(-b/6.0*(r*r)))*((1-r*r)**(b/3))

  Case (WLC)
     expphi = Exp( -b/6.0 * ( 2*r*r + 1/(1-r) -r -1 ));
  End Select
End Function expphi


Function  lam1_th(hs, N)
  Use bspglocons
  Implicit None
  Real (FPREC)  lam1_th
  Real (FPREC) ,  Intent (in) :: hs
  Integer,  Intent (in) :: N


  Real (FPREC) s,b,aj


  b   = 1 - 1.66*hs**0.78;
  s   =   - 1.40*hs**0.78;

  aj = 4*Sin(Pi/2./N)*Sin(Pi/2./N);
  lam1_th = 2./( aj * b / N**s);

End Function lam1_th


!!! Time-stamp: <gsipceq.f90 14:12, 22 Sep 2005 by P Sunthar>

!________________________________________________
!   Bead Spring Configuration Space Utilities
!________________________________________________

!!! $Log: gsipceq.f90,v $
!!! Revision 1.4  2005/05/24 01:06:20  pxs565
!!! With variable NTOPB
!!!
!!! Revision 1.3  2005/05/02 11:24:40  pxs565
!!! Tuned to work with intel/lapack
!!!
!!! Revision 1.2  2005/03/16 00:49:09  pxs565
!!! Preaveraged RPY working
!!!
!!! Revision 1.1  2005/02/25 10:21:22  pxs565
!!! Initial revision
!!!


Module csputls 
!!! This module contains subroutines used to manipulate the 
!!! configuration variable, R and its derivatives
  Use bspglocons

Contains

  Subroutine b2bvector_sym_up(N,R,b2b)
    Implicit None
    ! Calling arguments
    Integer, Intent(in)  :: N
    Real (FPREC) , Intent(in)  :: R(:,:)
    Real (FPREC) , Intent(out) :: b2b(:,:,:)

!!!_________________________________________________________|
!!! This routine returns the vector one bead to another bead
!!! as a symmetric matrix but only the upper elements filled
!!!_________________________________________________________|

    Integer mu,nu

    ! according to /usr/bin/prof, this takes max time. so commenting out
    ! since only upper diagonal is anyway reqd
    ! b2b = 0

    ! note only upper half (nu > mu) of the matrix is filled
    ! and the convention is R_12 = R_2 - R_1
    ! where R_12 is the vector from 1 to 2
    Do nu = 2,N
       Do mu = 1,nu-1
          b2b(:,mu,nu) = R(:,nu) - R(:,mu)
       End Do
    End Do
  End Subroutine b2bvector_sym_up

  Subroutine modr_sym_up(N,b2bvec,deltaR) 
    Implicit None
    ! calling arguments
    Integer, Intent (in)  :: N
    Real (FPREC) , Intent (in)  :: b2bvec(:,:,:)
    Real (FPREC) , Intent (out) :: deltaR(:,:) 
!!!_________________________________________________________|
!!! This subroutine returns the magnitude of each of the bead 
!!! to bead vector
!!!_________________________________________________________|

    Integer mu,nu,i
    Real (FPREC)  r12(Ndim),modr


    deltaR = 0

    ! note that we cant use a forall since dot_product is a
    ! transformational function
    Do nu = 2,N
       Do mu = 1,nu-1
          r12 = b2bvec(:,mu,nu)
          modr = 0.0
          Do i=1,Ndim
             modr = modr + r12(i) * r12(i)
          End Do
          If (modr < 1.0e-12) modr = 1.0e-12   !so that initial dr is != 0
          deltaR(mu,nu) = Sqrt(modr)
       End Do
    End Do
  End Subroutine modr_sym_up


  Subroutine tensor_prod_of_vec(alpha,vec,tensor)
    Implicit None
    Real (FPREC) , Intent (in) :: alpha
    Real (FPREC) , Intent (in) :: vec(:)
    Real (FPREC) , Intent (out) ::  tensor(:,:)
    Integer i,j
    Real (FPREC)  temp

    Do j = 1,Ndim
       Do i = 1,j
          temp = alpha * vec(i) * vec(j)
          tensor(i,j) = temp
          tensor(j,i) = temp
       End Do
    End Do
  End Subroutine tensor_prod_of_vec
End Module csputls



Module BSpModel
!!! This module contains the constants and functions which are specific
!!! to the model of spring and EV.
  Use bspglocons
  Implicit None

  ! These are paramters that will be used  by the following procedures
  Real (FPREC)  sqrtb

  ! paramters for EV
  Real (FPREC)  zsbyds5,pt5bydssq

  ! preaveraged HI
  !Real (FPREC), save, dimension(Ndim,NBeads,Ndim,NBeads) ::  Bpmat,Dpmat
  Real (FPREC), save, Allocatable, Dimension(:,:,:,:) ::  Bpmat,Dpmat


Contains 

  Subroutine force_sans_hookean(sptype,sr,ff)
    Integer, Intent (in) :: sptype 
    Real (FPREC) , Intent (in) :: sr
    Real (DOBL) , Intent (out) :: ff

    Real (DOBL) r

    r = sr 
    ! quick way to do double prec calculation.  otherwise, the complete 
    ! code needs to be Dprec to incorporate accuracy at every level

    ! compute the force factor, ff
    Select Case (sptype)
    Case (HOOK) ! Hookean
       ff = 1.0

    Case (FENE) ! Warner Spring
       ff = 1./(1 - r*r)

    Case (ILC)  ! Inverse Langevin, Pade approx 
       ff = (3-r*r)/(1 - r*r)/3.

    Case (WLC) ! Worm-like Chain, Marko-Siggia interpolation
       ff = 1./(6*r) * ( 4*r+ 1./(1-r)/(1-r) - 1)

    End Select
  End Subroutine force_sans_hookean

  Subroutine Spring_force(sptype,N,b2bvec,dr,Fs)
    Integer, Intent (in) :: sptype
    Integer, Intent (in)  :: N
    Real (FPREC) , Intent (in)  :: b2bvec(:,:,:)! Ndim x Nbeads x Nbeads
    Real (FPREC) , Intent (in)  :: dr(:,:)      ! Nbeads x Nbeads
    Real (FPREC) , Intent (out) :: Fs(:,:)      ! Ndim x Nbeads 
!!! Purpose:
    ! computes the net spring force on a bead due to 
    ! connector springs of neighbouring beads

    Integer nu
    Real (FPREC)  r ! = |Q| / sqrtb
    Real (DOBL)  Fcnu(Ndim), Fcnu1(Ndim), nonhook

    Fs  = 0.0

    Fcnu1 = 0.0 ! Fc of nu-1

    Do nu = 1,N-1 ! over each of the connector

       r = dr(nu,nu+1)/sqrtb ! only upper diagonal
       If (Abs(r - 1.0) < MYEPS) r = 1 - MYEPS

       Call force_sans_hookean(sptype,r,nonhook)

       ! connector force from 1 to 2 (nu to nu+1)
       Fcnu = b2bvec(:,nu,nu+1) * nonhook

       ! total spring force on a bead
       Fs(:,nu) = Fcnu - Fcnu1
       Fcnu1 = Fcnu ! save for the next bead
    End Do

    ! and the last bead
    Fs(:,N) = - Fcnu1

  End Subroutine Spring_force


  Subroutine solve_implicit_r(sptype,dtby4,gama,r,ff)
    Integer, Intent (in) :: sptype
    Real (FPREC) , Intent (in) :: dtby4, gama
    Real (FPREC) , Intent (out) :: r
    Real (DOBL) , Intent (out) :: ff

    !! Purpose
    ! solves the equation r ( 1 + dt/4 ff ) == gama
    ! where ff depends on spring force law F = H Q ff,
    ! for r and returns, r and ff(r)

    Real (FPREC)  coeff(4),denom

    ! set up the polynomial equation (cubic), and the guess for
    ! gama >> 1, obtained by the asymptotic behaviour

    coeff(4) = 1.0

    Select Case (sptype)
    Case (HOOK) ! Hookean
       r = gama/(1+dtby4)
       ff = 1
       Return

    Case (FENE) ! FENE
       coeff(1) = gama
       coeff(2) = -(1. + dtby4)
       coeff(3) = -gama
       r = 1 - dtby4/2/gama 

    Case (ILC) ! ILC
       denom = 3. + dtby4
       coeff(1) = 3. * gama / denom
       coeff(2) = -(3. + 3.* dtby4) / denom
       coeff(3) = -3. * gama / denom
       r = 1 - dtby4/3/gama

    Case (WLC) ! WLC
       denom = 2*(1.5 + dtby4)
       coeff(1) = -3 *gama/denom
       coeff(2) =  3 * (1. + dtby4 + 2.*gama) / denom
       coeff(3) = -1.5 * (4. + 3.*dtby4 + 2.*gama) / denom
       r = 1 - Sqrt(dtby4/6/gama)

    End Select

    ! the Hookean guess is same for all the force laws
    If (gama < 1.0) Then 
       r = gama/(1.+dtby4)  
    End If

    ! all the common forces laws yeild a cubic in the implicit form
    ! polish the guessed root by a Newt-Raph for polynomials

    ! for WLC, the newt raph, does not converge to the correct root,
    ! for gama > 100 , therefore simply ignore polishing
    If (sptype == WLC .And. gama > 100) Then
       r = r
    Else 
       Call polish_poly_root(coeff,r,1e-6)
    End If


    Call force_sans_hookean(sptype,r,ff)

  End Subroutine solve_implicit_r


  Subroutine Excluded_volume_force(N,b2bvec,dr,Fev,chbjerrium,&
  ldebye,zstar,rcstar,dstar,cutoffL)
    Use bspglocons

!!! This routine returns excluded volume force
    Integer, Intent (in) :: N
    Real (FPREC) , Intent (in)  :: b2bvec(:,:,:)! Ndim x Nbeads x Nbeads
    Real (FPREC) , Intent (in)  :: dr(:,:)      ! Nbeads x Nbeads
    Real (FPREC) , Intent (out) :: Fev(:,:)     ! Ndim x Nbeads 
! EV parameters (non-dim)
       Real (FPREC), Intent(in) :: zstar, rcstar, dstar,cutoffL
! Electrostatic parameters (non-dim)
       Real (FPREC), Intent(in) :: chbjerrium, ldebye


    Integer nu,mu
    Real (FPREC)  Fpair12(Ndim),pi2,lowerlimits

    Fev  = 0.0
    pi2  = PI*2.0/rcstar
    !lowerlimits=rcstar*.7127
    lowerlimits=0.7

       ! caution donot use forall here, i means F12 will b evaluated for 
       ! all the values of mu,nu and then assigned, which is not what we 
       ! want.  The rule to use forall statements is that lhs must be 
       ! unique for-all the values of loop index, and RHS must not depend 
       ! on any other index of lhs, other than the lhs itself
       ! ie, we can have  foo(mu) = foo(mu) + bar(nu)
       ! but not  foo(mu) = foo(mu+3) + bar(nu)
!     write(*,*)'zsbyds5,pt5bydssq'
!     write(*,*)zsbyds5,pt5bydssq
       Do nu = 2,N
          Do mu  = 1,nu-1 
             ! force between the pair mu,nu, since the EV force is
             ! repulsive, we take it to be in the opposite direction
             ! of the connector vector mu -> nu
             ! convention followed: force is positive for attraction
         Fpair12=0.d0
!    If (zsbyds5 .NE. 0.0) Then
!             Fpair12 = - b2bvec(:,mu,nu) & 
!                  * zsbyds5*Exp(-dr(mu,nu)*dr(mu,nu) * pt5bydssq) 
     if(abs(dr(mu,nu)).ge.rcstar)then
          Fpair12=0.d0
           elseif(abs(dr(mu,nu)).le.lowerlimits)then
        Fpair12 = - 4.0*b2bvec(:,mu,nu)*zstar/(lowerlimits*lowerlimits )*&
                 (12.0*(1.0/lowerlimits)**12-(6.0*(1.0/lowerlimits)**6))
         else

      Fpair12 = -4.0*b2bvec(:,mu,nu)*zstar/(dr(mu,nu)*dr(mu,nu))*&
              ((12.0*(1.0/dr(mu,nu))**12)-(6.0*(1.0/dr(mu,nu))**6))

    ! if(abs(dr(mu,nu)).ge.rcstar)then
    !     Fpair12=0.d0
    !    elseif(abs(dr(mu,nu)).le.lowerlimits)then
    !Fpair12 = - b2bvec(:,mu,nu)*zstar/(lowerlimits*lowerlimits )*&
    !          ((rcstar/lowerlimits)**12-(rcstar/lowerlimits)**6 &
    !          +dstar*lowerlimits*sin(pi2*lowerlimits))
     !
     !else
    !Fpair12 = - b2bvec(:,mu,nu)*zstar/(dr(mu,nu)*dr(mu,nu))*&
    !          ((rcstar/dr(mu,nu))**12-(rcstar/dr(mu,nu))**6 &
    !          +dstar*dr(mu,nu)*sin(pi2*dr(mu,nu)))
    !    endif
    End If
! Electric potential
        if(abs(dr(mu,nu)).le.lowerlimits)then
   Fpair12 = Fpair12 - b2bvec(:,mu,nu)*Exp(-lowerlimits/ldebye)*&
chbjerrium/(lowerlimits*lowerlimits)*(1.0/lowerlimits+1.0/ &
ldebye)
     else
   Fpair12 = Fpair12 - b2bvec(:,mu,nu)*Exp(-dr(mu,nu)/ldebye)*&
chbjerrium/(dr(mu,nu)*dr(mu,nu))*(1.0/dr(mu,nu)+1.0/ &
ldebye)


     endif
!      write(*,*) "Fpair1", Fpair12
        Fev(:,mu) = Fev(:,mu) + Fpair12
 !     write(*,*) "Fpair3",mu, Fev(:,mu)
             Fev(:,nu) = Fev(:,nu) - Fpair12
  !    write(*,*) "Fpair4",nu, Fev(:,nu)
          End Do
   !   write(*,*) "Fpair5",mu, Fev(:,mu)
    !  write(*,*) "Fpair6",nu, Fev(:,nu)
       End Do
     ! write(*,*) "Fpair7",mu,nu, Fev(:,:)
  End Subroutine Excluded_volume_force

  Subroutine preavDnB(N,hs)

    Integer , Intent (in) :: N
    Real (FPREC), Intent (in) :: hs

    Integer mu,nu,i,j,k
    Integer lwork, info
    Real (FPREC)  x, erfx, Hbar(N,N), Hp(N*(N+1)/2)  &
         , Evec(N,N), Eval(N) &
         , Bpreav(N,N)

    Real (FPREC), Allocatable :: work(:)


    Real (FPREC), save :: hsave=0,nsave=0

    if (hsave ==  hs .and. nsave == N) Then
       Return  ! to avoid recalculation of Bpmat and Dpmat
    end if

    If (.Not. Allocated(Bpmat) ) Then
       Allocate(Bpmat(Ndim,N,Ndim,N),  &
            Dpmat(Ndim,N,Ndim,n))
    Else if (nsave /= N) Then
       Deallocate(Bpmat,Dpmat)
       Allocate(Bpmat(Ndim,N,Ndim,N),  &
            Dpmat(Ndim,N,Ndim,n))
    End If

    hsave = hs
    nsave = N

    !-- Calculate the preaveraged diffusion tensor (only upper diagonal)
    If (hs == 0)  Then
       Write (*,*) 'Hstar = 0, called in preavDnB, check call sequence'
       Return
    end If

    Hbar = 0

    Do mu = 1,N
       Do nu = mu+1,N
          x = sqrt(2*Pi/(nu-mu)) * hs

          Select Case (FPREC)
          Case(SNGL)
          erfx = erf(x) ! plain fortran (calling libc ?)
             !call vserf(1,x,erfx) ! Intel fortran
             !erfx = erf(x) ! IBM XLF
          Case(DOBL)
          erfx = erf(x) ! plain fortran (calling libc ?)
             !call vderf(1,x,erfx) ! Intel fortran
             !erfx = derf(x) ! IBM XLF
             !--
             !  Note that with IBM, this line needs to be commented out
             !  for single precision, the compiler does not ignore the
             !  derf(real x) call
             !--
          end Select

          !-- For RPY tensor (Fixman, JCP 78, 1594)
          Hbar(mu,nu) =  erfx - (1-exp(-x*x))/sqrt(Pi)/x
       end Do

       Hbar(mu,mu) = 1
    end Do


    !-- eigensystem of Hbar

    !!*** IBM (ESSL) only ***!!
!!$    !-- 
!!$    !  ESSL: convert to upper packed storage mode, see example in
!!$    !   http://publib.boulder.ibm.com/clresctr/docs/essl/essl_42/
!!$    !   200410/am501402/am501402186.html
!!$    !--
    k=0
    Do j=1,N
       Do i=1,j
          k = k+1
          Hp(k) = Hbar(i,j)
       end Do
    end Do

    lwork = 0 ! automatic allocation
    Select Case (FPREC)
    Case(SNGL)
       call sspev(21,Hp,Eval,Evec,N,N,work,lwork)
    Case(DOBL)
       call dspev(21,Hp,Eval,Evec,N,N,work,lwork)
    end Select

    !!*** LAPACK  ***!!

!    
Comments: