subroutine WriteData( InputConf, Nham, Nljgrs, Niongrs, & Nsp, MaxSp, MaxSteps, MaxBeads, MaxInt, MaxReal, & BETA, ZETA, NAMElj, MASSlj, EPS, SIG, CP, ALP, KAPPA, & LAMDA, NAMEion, MASSion, CHARGE, NAMEsp, NSTEPS, & LENLJ, LENION, BEADTYPE, GROUPTYPE, NTRIALS, & STEPSTART, STEPLENGTH, METHOD, INTPARAM, REALPARAM, & ERSTEPS, ERSTART, EREND, EESTEPS, MaxEE, BEADDAMP, & FromDisk, Nmol, Volume, nmoves, PROB_MOVE, PROB_DR, & PROB_SP_CD, PROB_SP_RG, NinCycle, Nstorage, Nadjust, & Ncollect, Uwidth, Seed, ResultsFile ) implicit none ! InputConf has the information for continuing a previous run. character*30, intent(in) :: InputConf ! 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) :: 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. ! mu is the chemical potential term real, dimension(Nham), intent(in) :: BETA real, dimension(Nsp,Nham), intent(in) :: ZETA ! 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 ! 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. ! 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 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 integer, dimension(MaxSp) :: ERSTEPS integer, dimension(MaxSteps,MaxSp) :: ERSTART integer, dimension(MaxSteps,MaxSp) :: EREND integer, dimension(MaxSp) :: EESTEPS integer :: MaxEE real, dimension(MaxBeads,MaxEE,MaxSp) :: BEADDAMP ! 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), intent(in) :: Nmol ! BoxSize contains the length of the simulation boxes. real, intent(in) :: Volume ! 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_CD 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. integer, intent(in) :: nmoves real, dimension(nmoves), intent(in) :: PROB_MOVE real, dimension(2), intent(in) :: PROB_DR real, dimension(MaxSp), intent(in) :: PROB_SP_CD 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. ! Ncollect is the number of MC steps before collecting histogram data. integer, intent(in) :: NinCycle integer, intent(in) :: Nstorage integer, intent(in) :: Nadjust integer, intent(in) :: Ncollect ! Uwidth is the width of an energy bin. real, intent(in) :: Uwidth ! Seed is the current random number generator seed value. integer, intent(in) :: Seed ! ResultsFile is the file name to write the data to character*30, intent(in) :: ResultsFile ! Local stuff integer :: h, i, j, k, m, n integer, parameter :: unit = 60 real :: tmp real, parameter :: Pi = 3.14159265359 real, parameter :: Nav = 6.02214e23 real, external :: zbrent3, funk3 open(unit, file = ResultsFile, position = 'append' ) write(unit,"(1x,A)") 'EHSGCMC Simulation' write(unit,12 ) 'Temperature = ', 1.0/BETA, ' K ' 12 format( 1x, A, F10.3, A ) do i = 1, Nsp write(unit,13) 'ln(Zeta)', i, ' = ', log(ZETA(i,1:Nham)) 13 format( 1x, A, I2, A, F10.3 ) end do write(unit,"(1x,A, F10.3, A)" ) 'Volume = ', Volume, ' A^3 ' write(unit,*) write(unit,"(1x,A)") 'Exp-6 / LJ Groups: ' write(unit,"(1x,A)") 'Name Mass Hamil. Eps/k_B (K) Sigma (A) C Alpha ' do h = 1, Nham do i = 1, Nljgrs if( CP(i,i,h) > 0.0 ) then tmp = SIG(i,i,h) * zbrent3( funk3, 0.5, 1.0, 1.0e-7, ALP(i,i,h), CP(i,i,h) ) else tmp = SIG(i,i,h) end if write(unit,"(1x,A,T13,f8.4,T24,I2,T31,f7.3,T44,f7.4,T54,f7.4,T65,f5.2)") NAMElj(i), MASSlj(i), h, & EPS(i,i,h), tmp, CP(i,i,h), ALP(i,i,h) end do end do do h = 1, Nham write(unit, "(1x,A,I2)") 'Cross Parameters for hamiltonian ', h write(unit, "(T21,A,T51,A)") 'Kappa_ij', 'Lamda_ij' write(unit, 10) ' Group i\j', (i, i=1,Nljgrs), (i, i=1,Nljgrs) 10 format( A,T14,I8,T44,I8 ) do i = 1, Nljgrs write(unit, 20) i, (KAPPA(i,j,h),j=1,Nljgrs) , (LAMDA(i,j,h),j=1,Nljgrs) 20 format( I7,T16,F8.4,T46,F8.4 ) end do end do if( Niongrs > 0 ) then write(unit,*) write(unit,"(1x,A)") 'Ionic Groups: ' write(unit,"(1x,A)") 'Name Mass Hamil. Charge (e) ' do h = 1, Nham do i = 1, Niongrs write(unit,"(1x,A,T13,G12.5,T27,I2,T38,SP,F8.5,SS)") NAMEion(i), MASSion(i), h, & 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 = ', Nmol(i) 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, 31) & INTPARAM(3:2+m,j,i), & 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 31 format( T30, I3, T45, <1+2*m>G11.5 ) 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, 32) & INTPARAM(5:4+m+2*n,j,i) if( INTPARAM(3,j,i) == 1 ) then write(unit, 33) & REALPARAM(1:1+2*m+6*n,j,i) else if( INTPARAM(3,j,i) == 2 ) then write(unit, 34) & 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 32 format( T30, I3 ) 33 format( T45, <1+2*m+6*n>G11.5 ) 34 format( T45, <1+2*m+4*n>G11.5 ) 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, 35) & INTPARAM(3:2+m,j,i), & 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 35 format( T30, I3, T45, <2+2*m>G11.5 ) 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, 36) & INTPARAM(5:4+m+2*n,j,i) if( INTPARAM(3,j,i) == 1 ) then write(unit, 37) & REALPARAM(1:2+2*m+6*n,j,i) else if( INTPARAM(3,j,i) == 2 ) then write(unit, 38) & 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 36 format( T30, I3 ) 37 format( T45, <2+2*m+6*n>G11.5 ) 38 format( T45, <2+2*m+4*n>G11.5 ) 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)") ' Configurational Bias Steps ' 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 write(unit,"(A)") ' Early Rejection Steps' write(unit,"(A)") ' ER Step Start CCB Step End CCB Step ' do j = 1, ERSTEPS(i) write(unit, "(1x,I3,T18,I3,T36,I3)") j, ERSTART(j,i), EREND(j,i) end do write(unit,"(A)") ' Expanded Ensemble Damping Factors' write(unit,"(A,I2)") ' Bead Alpha 1 to ', EESTEPS(i) do j = 1, LENLJ(i) + LENION(i) write(unit,"(1x,I3,T12, 10(f7.4,2x))") j, BEADDAMP(j,1:EESTEPS(i),i) 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 write(unit,"(1x,A,T48,A,I10)") 'Initial Seed', '= ', Seed write(unit,"(1x,A,T48,A,4F7.2)") 'Percentage of Disp/Rot, Creat/Destruc, Regrow', & '= ', PROB_MOVE(1) * 100.0, & ( PROB_MOVE(2) - PROB_MOVE(1) ) * 100.0, & ( PROB_MOVE(3) - PROB_MOVE(2) ) * 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 Creat/Destruc', '= ' do i = 1, Nsp write(unit,"(F7.2)", Advance='No') ( PROB_SP_CD(i) - sum( PROB_SP_CD(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,"(1x,A,T48,A,I8)") 'MC Steps per Histogram data Collection', '= ', & Ncollect write(unit,"(1x,A,T48,A,F10.4)") 'Width of energy histogram bin (K)', '= ', & Uwidth write(unit,*) close(unit) return end subroutine WriteData