subroutine WriteData( Ensemble, InputConf, Ntemp, Nham, Nljgrs, Niongrs, & Nsp, MaxSp, MaxSteps, MaxBeads, MaxInt, MaxReal, & BETA, Pressure, NAMElj, MASSlj, EPS, SIG, CP, ALP, & KAPPA, LAMDA, GAMMA, NAMEion, MASSion, CHARGE, NAMEsp, & NSTEPS, LENLJ, LENION, BEADTYPE, GROUPTYPE, MASSsp, & NTRIALS, STEPSTART, STEPLENGTH, METHOD, INTPARAM, & REALPARAM, FromDisk, Nmol, BoxSize, PROB_MOVE, PROB_DR, & PROB_SP_TR, PROB_SP_RG, NinCycle, Nstorage, Nadjust, & Seed, unit ) implicit none ! Ensemble is the type of Gibbs-ensemble being used. ! InputConf has the information for continuing a previous run. character*5, intent(in) :: Ensemble character*30, intent(in) :: InputConf ! Ntemp is the number of temperatures. ! Nham is the number of hamiltonians. ! Nljgrs is the number of Lennard-Jones groups in the simulation. ! Niongrs is the number of ionic groups in the simulation. ! Nsp is the number of species in the simulation. integer, intent(in) :: Ntemp integer, intent(in) :: Nham integer, intent(in) :: Nljgrs integer, intent(in) :: Niongrs integer, intent(in) :: Nsp ! MaxSp is the maximum number of species in the simulation. ! MaxSteps is the maximum number of CB steps for a molecule. ! MaxBeads is the maximum number of LJ and ionic beads in a molecule. ! MaxInt is the maximum number of integer parameters for a CB method. ! MaxReal is the maximum number of real parameters for a CB method. integer, intent(in) :: MaxSp integer, intent(in) :: MaxSteps integer, intent(in) :: MaxBeads integer, intent(in) :: MaxInt integer, intent(in) :: MaxReal ! BETA is the reciprical temperature. ! Pressure is the pressure for a constant pressure run. real, dimension(Ntemp,1), intent(in) :: BETA real, intent(in) :: Pressure ! NAMElj containes the names of the LJ groups. ! MASSlj contains the mass of the LJ groups. ! EPS contains the eps_ij parameters for each hamiltonian. ! SIG contains the sigma_ij parameters for each hamiltonian. ! CP contains the C_ij parameters for each hamiltonian. ! ALP contains the alpha_ij parameters for each hamiltonian. ! KAPPA contains the correction to the epsilom_ij cross parameters. ! LAMDA contains the correction to the sigma_ij cross parameters. character*15, dimension(Nljgrs), intent(in) :: NAMElj real, dimension(Nljgrs), intent(in) :: MASSlj real, dimension(Nljgrs,Nljgrs,Nham), intent(in) :: EPS real, dimension(Nljgrs,Nljgrs,Nham), intent(in) :: SIG real, dimension(Nljgrs,Nljgrs,Nham), intent(in) :: CP real, dimension(Nljgrs,Nljgrs,Nham), intent(in) :: ALP real, dimension(Nljgrs,Nljgrs,Nham), intent(in) :: KAPPA real, dimension(Nljgrs,Nljgrs,Nham), intent(in) :: LAMDA real, dimension(Nljgrs,Nljgrs,Nham), intent(in) :: GAMMA ! NAMEion containes the names of the ionic groups. ! MASSion contains the mass of the ionic groups. ! CHARGE contains the charge of a bead for each hamiltonian. character*15, dimension(Niongrs), intent(in) :: NAMEion real, dimension(Niongrs), intent(in) :: MASSion real, dimension(Niongrs,Nham), intent(in) :: CHARGE ! NAMEsp contains the name of each species. ! NSTEPS is the number of CB steps needed to grow each species. ! LENLJ is the LJ length of each species. ! LENION is the ionic length of each species. ! BEADTYPE indicates whether a bead is LJ or ionic. ! GROUPTYPE indicates the group identity of each bead. ! MASSsp contains the mass of the species. ! NTRIALS is the number of trials for each CB step. ! STEPSTART is the starting bead for each CB step. ! STEPLENGTH is the bead length of each CB step. ! METHOD is the method used for each CB step. ! INTPARAM contains the integer parameters for each CB step. ! REALPARAM contains the real parameters for each CB step. character*15, dimension(Nsp), intent(in) :: NAMEsp integer, dimension(MaxSp), intent(in) :: NSTEPS integer, dimension(MaxSp), intent(in) :: LENLJ integer, dimension(MaxSp), intent(in) :: LENION character*5, dimension(MaxBeads,MaxSp), intent(in) :: BEADTYPE integer, dimension(MaxBeads,MaxSp), intent(in) :: GROUPTYPE real, dimension(MaxSp), intent(in) :: MASSsp integer, dimension(MaxSteps,MaxSp), intent(in) :: NTRIALS integer, dimension(MaxSteps,MaxSp), intent(in) :: STEPSTART integer, dimension(MaxSteps,MaxSp), intent(in) :: STEPLENGTH character*10, dimension(MaxBeads,MaxSp), intent(in) :: METHOD integer, dimension(MaxInt,MaxBeads,MaxSp), intent(in) :: INTPARAM real, dimension(MaxReal,MaxBeads,MaxSp) :: REALPARAM ! FromDisk indicates whether or not to read the initial configuration from disk. logical, intent(in) :: FromDisk ! Nmol is the total and per species number of molecules in the simulation. integer, dimension(0:MaxSp, 0:2), intent(in) :: Nmol ! BoxSize contains the length of the simulation boxes. real, dimension(2), intent(in) :: BoxSize ! PROB_MOVE is the accumulative probability of performing a ! Displacement / Rotation, Volume change, transfer, or regrowth. ! PROB_DR is the accumulative probability of performing a Displacement or Rotation. ! PROB_SP_TR is the accumulative probability of selecting a given species for transfer. ! PROB_SP_RG is the accumulative probability of selecting a given species for regrowth. real, dimension(4), intent(in) :: PROB_MOVE real, dimension(2), intent(in) :: PROB_DR real, dimension(MaxSp), intent(in) :: PROB_SP_TR real, dimension(MaxSp), intent(in) :: PROB_SP_RG ! NinCycle is the number of MC steps per cycle. ! Nstorage is the number of MC steps before storing the current configuration. ! Nadjust is the number of MC steps before adjusting the maximum displacement, etc. integer, intent(in) :: NinCycle integer, intent(in) :: Nstorage integer, intent(in) :: Nadjust ! Seed is the current random number generator seed value. integer, intent(in) :: Seed ! unit is the disk number to wright output to integer, intent(in) :: unit ! Local stuff integer :: h, i, j, k, m, n real :: bx1, bx2, tmp real, parameter :: Pi = 3.14159265359 real, parameter :: Nav = 6.02214e23 real, external :: zbrent3, funk3 write(unit,"(1x,A,A)") trim(Ensemble), '-Gibbs Ensemble Simulation' if( Ensemble == 'NPT' ) then write(unit,"(1x,A, F10.3, A)") 'Pressure = ', Pressure, ' bar' end if write(unit,"(1x,A, F10.3, A)" ) 'Temperature = ', 1.0/BETA(1,1), ' K ' write(unit,*) if( CP(1,1,1) > 0.0 ) then write(unit,"(1x,A)") 'Exp-6 Groups: ' write(unit,"(1x,A)") 'Name Mass Eps/k_B (K) Sigma (A) C Alpha' do h = 1, Nham do i = 1, Nljgrs tmp = SIG(i,i,h) * zbrent3( funk3, 0.5, 1.0, 1.0e-7, ALP(i,i,h), CP(i,i,h) ) write(unit,"(1x,A,T13,f9.5,T24,f8.3,T38,f7.4,T52,f7.4,T66,f6.2)") NAMElj(i), MASSlj(i), & EPS(i,i,h), tmp, CP(i,i,h), ALP(i,i,h) end do end do else write(unit,"(1x,A)") 'Lennard-Jones Groups: ' write(unit,"(1x,A)") 'Name Mass Eps/k_B (K) Sigma (A) ' do h = 1, Nham do i = 1, Nljgrs write(unit,"(1x,A,T13,f9.5,T24,f8.3,T38,f7.4)") NAMElj(i), MASSlj(i), & EPS(i,i,h), SIG(i,i,h) end do end do end if do h = 1, Nham write(unit, "(1x,A)") 'Cross Parameters ' if( Nljgrs < 4 ) then write(unit, "(T21,A,T46,A,T71,A)") 'Kappa_ij', 'Lamda_ij', 'Gamma_ij' write(unit, 10) ' Group i\j', (i, i=1,Nljgrs), (i, i=1,Nljgrs), (i, i=1,Nljgrs) 10 format( A, T14, I8, T39, I8, T64, I8 ) !10 format( A, T14, 4I8, T44, 4I8, T64, 4I8 ) do i = 1, Nljgrs write(unit, 20) i, (KAPPA(i,j,h),j=1,Nljgrs), (LAMDA(i,j,h),j=1,Nljgrs), & (GAMMA(i,j,h),j=1,Nljgrs) end do 20 format( I7, T16, F8.4, T41, F8.4, T66, F8.4 ) !20 format( I7, T16, 4F8.4, T46, 4F8.4, T64, 4I8 ) else write(unit, "(T21,A)") 'Kappa_ij' write(unit, "(A,T14,20I8)") ' Group i\j', (i, i=1,Nljgrs) do i = 1, Nljgrs write(unit, "(I7,T16,20F8.4)") i, (KAPPA(i,j,h),j=1,Nljgrs) end do write(unit, "(T21,A)") 'Lamda_ij' write(unit, "(A,T14,20I8)") ' Group i\j', (i, i=1,Nljgrs) do i = 1, Nljgrs write(unit, "(I7,T16,20F8.4)") i, (LAMDA(i,j,h),j=1,Nljgrs) end do write(unit, "(T21,A)") 'Gamma_ij' write(unit, "(A,T14,20I8)") ' Group i\j', (i, i=1,Nljgrs) do i = 1, Nljgrs write(unit, "(I7,T16,20F8.4)") i, (GAMMA(i,j,h),j=1,Nljgrs) end do end if end do if( Niongrs > 0 ) then write(unit,*) write(unit,"(1x,A)") 'Ionic Groups: ' write(unit,"(1x,A)") 'Name Mass Charge (e) ' do h = 1, Nham do i = 1, Niongrs write(unit,"(1x,A,T13,G12.5,T27,F8.5)") NAMEion(i), MASSion(i), & CHARGE(i,h) end do end do end if do i = 1, Nsp write(unit,*) write(unit,"(1x,A,A)") trim( NAMEsp(i) ), ' Data' write(unit,"(1x,A,I6,I6)") 'Starting number of molecules in 1, 2 = ', Nmol(i,1), Nmol(i,2) write(unit,"(1x,A,A)") 'Growth of ', NAMEsp(i) write(unit,"(1x,A)") 'Bead Type Group Method Int. Params Real Params' do j = 1, LENLJ(i) + LENION(i) select case ( METHOD(j,i) ) case( 'Random' ) write(unit,"(1x,I3,T9,A,T13,I3,T21,A)") & j, BEADTYPE(j,i), GROUPTYPE(j,i), METHOD(j,i) case( 'Sphere' ) write(unit,"(1x,I3,T9,A,T13,I3,T21,A,T30,I3,T45,G11.5)") & j, BEADTYPE(j,i), GROUPTYPE(j,i), METHOD(j,i), & INTPARAM(1,j,i), REALPARAM(1,j,i) case( 'Cone' ) write(unit,"(1x,I3,T9,A,T13,I3,T21,A,T30,3I3,T45,20(G11.5,1x))") & j, BEADTYPE(j,i), GROUPTYPE(j,i), METHOD(j,i), & INTPARAM(1:3,j,i), REALPARAM(1:2*INTPARAM(3,j,i),j,i) * 180.0/Pi, & REALPARAM(2*INTPARAM(3,j,i)+1:3*INTPARAM(3,j,i),j,i) case( 'ConeTor' ) if( INTPARAM(4,j,i) == 1 ) then write(unit,"(1x,I3,T9,A,T13,I3,T21,A,T30,4I3,T45,9G11.5)") & j, BEADTYPE(j,i), GROUPTYPE(j,i), METHOD(j,i), & INTPARAM(1:4,j,i), & REALPARAM(1,j,i) * 180.0/Pi, REALPARAM(2:8,j,i) else if( INTPARAM(4,j,i) == 2 ) then write(unit,"(1x,I3,T9,A,T13,I3,T21,A,T30,4I3,T45,7G11.5)") & j, BEADTYPE(j,i), GROUPTYPE(j,i), METHOD(j,i), & INTPARAM(1:4,j,i), & REALPARAM(1,j,i) * 180.0/Pi, REALPARAM(2:6,j,i) end if case( 'FxBend' ) write(unit,"(1x,I3,T9,A,T13,I3,T21,A,T30,2I3)") & j, BEADTYPE(j,i), GROUPTYPE(j,i), METHOD(j,i), & INTPARAM(1:2,j,i) m = INTPARAM(1,j,i) do k = 1, m REALPARAM( 1+2*k, j, i ) = REALPARAM( 1+2*k, j, i ) * 180.0 / Pi end do write(unit, "(T30, 5I3)") INTPARAM(3:2+m,j,i) write(unit, "(T45, 10G11.5)") REALPARAM(1:1+2*m,j,i) do k = 1, m REALPARAM( 1+2*k, j, i ) = REALPARAM( 1+2*k, j, i ) * Pi / 180.0 end do case( 'FxBendTor' ) write(unit,"(1x,I3,T9,A,T13,I3,T21,A,T30,4I3)") & j, BEADTYPE(j,i), GROUPTYPE(j,i), METHOD(j,i), & INTPARAM(1:4,j,i) m = INTPARAM(1,j,i) n = INTPARAM(2,j,i) do k = 1, m REALPARAM( 1+2*k, j, i ) = REALPARAM( 1+2*k, j, i ) * 180.0 / Pi end do write(unit, "(T30, 20I3)") INTPARAM(5:4+m+2*n,j,i) if( INTPARAM(3,j,i) == 1 ) then write(unit, "(T45, 50G11.5)") REALPARAM(1:1+2*m+6*n,j,i) else if( INTPARAM(3,j,i) == 2 ) then write(unit, "(T45, 50G11.5)") REALPARAM(1:1+2*m+4*n,j,i) end if do k = 1, m REALPARAM( 1+2*k, j, i ) = REALPARAM( 1+2*k, j, i ) * Pi / 180.0 end do case( 'Stretch' ) write(unit,"(1x,I3,T9,A,T13,I3,T21,A,T30,I3,T45,2G11.5)") & j, BEADTYPE(j,i), GROUPTYPE(j,i), METHOD(j,i), & INTPARAM(1,j,i), REALPARAM(1:2,j,i) case( 'StBend' ) write(unit,"(1x,I3,T9,A,T13,I3,T21,A,T30,2I3)") & j, BEADTYPE(j,i), GROUPTYPE(j,i), METHOD(j,i), & INTPARAM(1:2,j,i) m = INTPARAM(1,j,i) do k = 1, m REALPARAM( 2+2*k, j, i ) = REALPARAM( 2+2*k, j, i ) * 180.0 / Pi end do write(unit, "(T30, 5I3)") INTPARAM(3:2+m,j,i) write(unit, "(T45, 10G11.5)") REALPARAM(1:2+2*m,j,i) do k = 1, m REALPARAM( 2+2*k, j, i ) = REALPARAM( 2+2*k, j, i ) * Pi / 180.0 end do case( 'StBendTor' ) write(unit,"(1x,I3,T9,A,T13,I3,T21,A,T30,4I3)") & j, BEADTYPE(j,i), GROUPTYPE(j,i), METHOD(j,i), & INTPARAM(1:4,j,i) m = INTPARAM(1,j,i) n = INTPARAM(2,j,i) do k = 1, m REALPARAM( 2+2*k, j, i ) = REALPARAM( 2+2*k, j, i ) * 180.0 / Pi end do write(unit, "(T30, 15I3)") INTPARAM(5:4+m+2*n,j,i) if( INTPARAM(3,j,i) == 1 ) then write(unit, "(T45, 50G11.5)") & REALPARAM(1:2+2*m+6*n,j,i) else if( INTPARAM(3,j,i) == 2 ) then write(unit, "(T45, 50G11.5)") & REALPARAM(1:2+2*m+4*n,j,i) end if do k = 1, m REALPARAM( 2+2*k, j, i ) = REALPARAM( 2+2*k, j, i ) * Pi / 180.0 end do case( 'Match' ) write(unit,"(1x,I3,T9,A,T13,I3,T21,A,T30,I3)") & j, BEADTYPE(j,i), GROUPTYPE(j,i), METHOD(j,i), & INTPARAM(1,j,i) case( 'Complete' ) write(unit,"(1x,I3,T9,A,T13,I3,T21,A)") & j, BEADTYPE(j,i), GROUPTYPE(j,i), METHOD(j,i) end select end do write(unit,"(A)") ' Step Attempts Start Bead End Bead' do j = 1, NSTEPS(i) write(unit, "(1x,I3,T12,I3,T26,I3,T40,I3)") j, NTRIALS(j,i), & STEPSTART(j,i), STEPSTART(j,i) + STEPLENGTH(j,i) - 1 end do end do write(unit,*) if( FromDisk ) then write(unit,"(1x,A,A)") 'The initial configuration was taken from file: ', & trim( InputConf ) else write(unit,"(1x,A)") 'The initial configuration was randomly generated ' end if bx1 = BoxSize(1) bx2 = BoxSize(2) if( bx1 > -1.0e-7 .AND. bx1 < 1.0e-7 ) bx1 = 1.0 if( bx2 > -1.0e-7 .AND. bx2 < 1.0e-7 ) bx2 = 1.0 write(unit,"(1x,A,T48,A,I10)") 'Initial Seed', '= ', Seed write(unit,"(1x,A,T48,A,2F7.2)") 'Initial Densities (kg/m^3)', '= ', & dot_product( Nmol(1:Nsp,1) , MASSsp(1:Nsp) ) * 1.0e27 / Nav / bx1**3 , & dot_product( Nmol(1:Nsp,2) , MASSsp(1:Nsp) ) * 1.0e27 / Nav / bx2**3 write(unit,"(1x,A,T48,A,4F7.2)") 'Percentage of Disp/Rot, Vol, Trans, Regrow', & '= ', PROB_MOVE(1) * 100.0, & ( PROB_MOVE(2) - PROB_MOVE(1) ) * 100.0, & ( PROB_MOVE(3) - PROB_MOVE(2) ) * 100.0, & ( PROB_MOVE(4) - PROB_MOVE(3) ) * 100.0 write(unit,"(1x,A,T48,A,2F7.2)") 'Percentage of Displacement, Rotations', & '= ', PROB_DR(1) * 100.0, & ( 1.0 - PROB_DR(1) ) * 100.0 if( Nsp > 1 ) then write(unit,"(1x,A,T48,A)", Advance='No') 'Percentage of Species Transfers', '= ' do i = 1, Nsp write(unit,"(F7.2)", Advance='No') ( PROB_SP_TR(i) - sum( PROB_SP_TR(1:i-1) ) ) * 100.0 end do write(unit,*) write(unit,"(1x,A,T48,A)", Advance='No') 'Percentage of Species Regrowths', '= ' do i = 1, Nsp write(unit,"(F7.2)", Advance='No') ( PROB_SP_RG(i) - sum( PROB_SP_RG(1:i-1) ) ) * 100.0 end do write(unit,*) end if write(unit,"(1x,A,T48,A,3(I8,2x))") 'MC Steps per Cycle, Storage, Adjustment', '= ', & NinCycle, Nstorage, Nadjust write(unit,*) return end subroutine WriteData