[table of contents] [master index] [comments] [modules] [programs] [variables] [types] [procedures]
PROGRAM
pscf_pd
PURPOSE
Main program for polymer self-consistent field theory (PSCF) for
spatially periodic (PD) microstructures. The program can treat
incompressible mixtures of any number of linear multiblock
copolymers, linear homopolymers, and small molecule solvents.
This program contains a main loop that reads an input script from
standard input. The loop, and the operation of the program, ends
when a line containing the word FINISH is encountered. The manual
contains detailed description of the format of this script, but a
brief description is also given here:
Input Script Format:
The first line of the input script must be of the "format i j", in
which i and j are major and minor version numbers for the input
script file format (e.g., "format 1 0" for v1.0). The rest of the
input script consists of a series of blocks of input data. Each
block begins with a blank line followed by a line containing a
capitalized operation flag (op_flag) string, such as 'CHEMISTRY',
'UNIT_CELL', 'DISCRETIZATION', 'ITERATION', etc. The remainder
of each block (if any) is a series of variable values. Each
variable is preceded by a comment line containing the name of
the variable, as used in the program. One and two-dimensional
array variables may be input in any of several formats (e.g., all
values on a single line, or on multiple lines). See the manual
for the format of such variables. Some operation flags, such
as 'ITERATE' and 'SWEEP', trigger an action, beyond reading in
variables. Reading of the script, and the program, ends when the
'FINISH' flag is encountered.
Within each block, input variables are read using the input()
interface of the io_mod module. This is an overloaded interface
that can be used to read in a scalar or array of integer, real,
character, or boolean data. The main program uses the comment
style 'A' (for 'Above') of io_mod, in which the value of each
input variable is preceded by a line containing the name of
the variable. Within each block, variables must appear in a
predetetermined format. The required format is documented in
the manual, and in example scripts provided with the program.
All data is read using the fortran read(5,*) format, so that
the data format and use of white space is flexible, as long
as each data value is of the expected data type.
SOURCE
program pscf_pd
use const_mod
use string_mod
use io_mod
use version_mod
use field_io_mod
use chemistry_mod, only : input_chemistry, rescale_vref, &
input_monomers, input_chains, input_solvents,&
input_composition, input_interaction, &
output_monomers, output_chains, output_solvents,&
output_composition, output_interaction, &
N_monomer, N_chain, N_block, N_solvent, chi
use unit_cell_mod, only : input_unit_cell, output_unit_cell, &
N_cell_param, cell_param, &
make_unit_cell, R_basis, G_basis
use group_mod, only : output_group
use grid_mod, only : ngrid, input_grid, allocate_grid, make_ksq
use basis_mod, only : N_wave, N_star, group, &
make_basis, output_waves, release_basis
use fft_mod, only : fft_plan, create_fft_plan, fft, ifft
use grid_basis_mod
use chain_mod
use scf_mod, only : density_startup, density
use iterate_mod, only : input_iterate_param, output_iterate_param, &
itr_algo, domain, &
iterate_NR_startup, iterate_NR, &
iterate_AM_startup, iterate_AM
use sweep_mod
use response_mod, only : response_startup, response_sweep
implicit none
! FFT variable
type(fft_plan) :: plan
! SCFT variables
real(long) :: omega(:,:) ! chemical potential field
! omega(monomer,basis function)
real(long) :: rho(:,:) ! monomer density field
! rho(monomer,basis function)
real(long) :: f_Helmholtz ! free energy, units kT / monomer
real(long) :: pressure ! pressure * V_monomer / kT
real(long) :: stress(:) ! d(f_Helmholtz)/d(cell_param)
allocatable :: omega, rho, stress
! Input parameters (all others declared in modules)
character(60) :: group_name ! name of crystal space group
real(long) :: chain_step ! contour length step size
real(long) :: vref_scale ! Factor to rescale reference volume
! Input and output file names
character(60) :: input_prefix ! prefix for input omega file:
! input_prefix//omega
character(60) :: output_prefix ! prefix for output files after iteration:
! output_prefix//out
! output_prefix//omega
! output_prefix//rho
! output_prefix//group
! output_prefix//waves
character(60) :: input_filename ! name of input field file
character(60) :: output_filename ! name of output field file
! Variable for field transformations
integer :: i1, i2, i3, alpha
integer, allocatable :: grid_size(:)
complex(long), allocatable :: k_grid(:,:,:,:)
real(long), allocatable :: r_grid(:,:,:,:)
real(long), allocatable :: omega_basis(:,:)
real :: ff, qR
real(long) :: rnodes ! number of grid points
character(25) :: fmt
! Variables for iteration (fixed chemistry)
integer :: extr_order ! extrapolation order = 1
integer :: itr ! iteration index
real(long) :: error ! error = max(residual)
logical :: converge ! true if converged
! character(10) :: algo
! Variables for sweep (sequence of parameters)
integer :: i, j ! step indices
real(long) :: s_max ! maximum value of variable s
real(long) :: s ! continuation variable
real(long) :: step ! actual step size
real(long) :: step_unit ! unit step size
! Operation selection string from input script
character(60) :: op_string ! Operation selection string
! Logical operation flags - set true as each operation requested
! Note: These are listed in normal sequence within input file
logical :: monomer_flag = .FALSE. ! monomer data read
logical :: chain_flag = .FALSE. ! chain data read
logical :: solvent_flag = .FALSE. ! solvent data read
logical :: composition_flag = .FALSE. ! composition data read
logical :: interaction_flag = .FALSE. ! interaction data read
logical :: unit_cell_flag = .FALSE. ! unit_cell made
logical :: discretize_flag = .FALSE. ! grid and ds made
logical :: prefix_flag = .FALSE. ! io file prefixes read
logical :: basis_flag = .FALSE. ! symmetrized basis made
logical :: omega_flag = .FALSE. ! initial omega exists
logical :: iterate_flag = .FALSE. ! 1st iteration requested
logical :: output_flag = .FALSE. ! deferred iterate output
logical :: sweep_flag = .FALSE. ! sweep requested
logical :: rho_flag = .FALSE. ! initial_rho_exist
! Timing variables
real(long) :: start_time, basis_time, scf_time
real(long) :: rpa_time
! File Unit numbers (parameters)
integer, parameter :: out_unit = 21 ! output summary
integer, parameter :: field_unit = 22 ! omega and rho fields
integer :: ierr ! error msg for file io
! File format version numbers
type(version_type) :: version ! input script format
!------------------------------------------------------------------
call cpu_time(start_time)
! Set defaults for parameter I/O - see io_mod
call set_echo(1) ! echo inputs to standard out
call set_com_style('A','A','A') ! comments on line above data
call set_com_use('R') ! replace comment in echoed output
call set_io_units(i=5,o=6) ! set standard in and out units
! Read file format version from input script, echo to stdout
call input_version(version, 5)
call output_version(version, 6)
! Main operation loop - infinite loop
op_loop : do
! Read operation string from stdin
read(5,*)
read(5,*) op_string
! Echo to stdout
write(6,*)
call output(trim(op_string),f="N",j="L",o=6)
! Print any deferred output from preceding ITERATE command
if (output_flag) then
if (trim(op_string) == "SWEEP") then
call output_summary(trim(output_prefix)//'0.')
call output_fields(trim(output_prefix)//'0.')
else
call output_summary(output_prefix)
call output_fields(output_prefix)
end if
output_flag = .FALSE.
endif
select case(trim(op_string))
case ("MONOMERS")
! Input N_monomer and kuhn array
! See chemistry_mod and users manual
! Preconditions
if (monomer_flag) then
write(6,*) "Error: MONOMERS can only be read once"
exit op_loop
end if
monomer_flag = .TRUE.
call input_monomers(5,'F')
case ("CHAINS")
! Input N_chains, block_monomer, block_length
! See chemistry_mod and users manual
! Check preconditions (must know N_monomer first)
if (chain_flag) then
write(6,*) "Error: CHAINS can only be read once"
exit op_loop
end if
if (.not.monomer_flag) then
write(6,*) "Error: Must read MONOMERS before CHAINS"
exit op_loop
end if
chain_flag = .TRUE.
call input_chains(5,'F')
case ("SOLVENTS")
! Input N_solvents, solvent_monomer, solvent_size
! See chemistry_mod and users manual
! Check preconditions (must know N_monomer first)
if (solvent_flag) then
write(6,*) "Error: SOLVENTS can only be read once"
exit op_loop
end if
if (.not.monomer_flag) then
write(6,*) "Error: Must read MONOMERS before SOLVENTS"
exit op_loop
end if
if (composition_flag) then
write(6,*) "Error: Cannot read SOLVENTS after COMPOSITION"
exit op_loop
end if
solvent_flag = .TRUE.
call input_solvents(5,'F')
case ("COMPOSITION")
! Input ensemble, phi_chain and phi_solvent, or mu_chain and mu_solvent
! See chemistry_mod and users manual
! Check preconditions (must know N_chain and N_solvent first)
if ( .not. chain_flag ) then
write(6,*) &
"Error: Must read CHAINS before COMPOSITION"
exit op_loop
end if
iterate_flag = .FALSE.
composition_flag = .TRUE.
call input_composition(5,'F')
case ("INTERACTION")
! Input Flory-Huggins chi interaction parameters.
! Input interaction_type and:
! chi (if interaction_type = 'chi'), or
! chi_a, chi_b, & temperature (if interaction_type = 'chi_T')
! See chemistry_mod and users manual
! Check preconditions (must know N_monomer first)
if (.not.monomer_flag) then
write(6,*) "Error: Must read MONOMERS before INTERACTION"
exit op_loop
end if
iterate_flag = .FALSE.
interaction_flag = .TRUE.
call input_interaction(5,'F')
case ('UNIT_CELL')
! Check preconditions
if (unit_cell_flag) then
write(6,*) "Error: UNIT_CELL can only be read once"
exit op_loop
end if
unit_cell_flag = .TRUE.
! Read unit cell parameters (see unit_cell_mod)
call input_unit_cell(5,'F')
allocate(grid_size(dim))
! Construct initial unit cell (see unit_cell_mod)
call make_unit_cell
case ('DISCRETIZATION')
! Read spatial and contour length discretization of PDE
! Check preconditions (Needs N_monomer and unit_cell parameters)
if (discretize_flag) then
write(6,*) "Error: DISCRETIZATION can only be read once"
exit op_loop
end if
if (.not.monomer_flag) then
write(6,*) "Error: Must read MONOMERS before DISCRETIZATION"
exit op_loop
else if ( .not. unit_cell_flag ) then
write(6,*) "Error: Must read UNIT_CELL before DISCRETIZATION"
exit op_loop
end if
discretize_flag = .TRUE.
! Input ngrid = number of points in FFT grid in each direction
call input_grid
call allocate_grid(N_monomer)
!call input(extr_order,'extr_order')
extr_order = 1
call input(chain_step,'chain_step')
case ('BASIS')
! Construct symmetry-adapated basis functions. See basis_mod
! Check preconditions (needs unit cell and grid)
if (basis_flag) then
write(6,*) "Error: BASIS can only be read once"
exit op_loop
end if
if (.not. unit_cell_flag) then
write(6,*) "Error: Must read UNIT_CELL before BASIS"
exit op_loop
else if (.not.discretize_flag) then
write(6,*) "Error: Must read DISCRETIZATION before BASIS"
exit op_loop
end if
basis_flag = .TRUE.
! Read name of space group used to construct basis functions.
! The string group_name can be a space group symbol, a space
! group number, or the name of a file containing the elements
! of the group. See space_groups in module space_group_mod.
call input(group_name,'group_name')
! Construct basis functions (see basis_mod)
call make_basis(R_basis,G_basis,group_name,ngrid,grid_flag=.TRUE.)
! Output N_star (# of symmetrized basis functions) to stdout
call output(N_star,'N_star',o=6)
! Allocate omega, rho, stress (internal routine)
call allocate_scf_arrays
case ('RESCALE')
! This command should be read immediately before iterate.
! Rescale monomer reference volume. Change values of kuhn,
! chi, block_length, and solvent_size parameter arrays, and
! the omega field, to obtain an equivalent set of parameters
! and fields.
! Check preconditions
if (.not.composition_flag) then
write(6,*) "Error: Must read COMPOSITION before RESCALE"
exit op_loop
else if (.not.interaction_flag) then
write(6,*) "Error: Must read INTERACTION before RESCALE"
exit op_loop
else if (.not.basis_flag) then
write(6,*) "Error: Must read BASIS before RESCALE"
exit op_loop
end if
iterate_flag = .FALSE.
omega_flag = .TRUE.
! Read name of input omega file
call input(input_filename, 'input_filename')
! Read scale factor: vref -> vref/vref_scale
call input(vref_scale, 'vref_scale')
! Read omega field
open(unit=field_unit,file=trim(input_filename), status='old', iostat=ierr)
if (ierr/=0) stop "Error while opening omega file"
call input_field(omega, field_unit)
close(field_unit)
! Rescale kuhn, chi, block_length, solvent_size
! See chemistry_mod
call rescale_vref(vref_scale)
! Rescale omega field
omega = omega/vref_scale
case ('ITERATE')
! Iterate to convergence for one set of parameters
! See iterate_mod
! Check preconditions
if (.not.composition_flag) then
write(6,*) "Error: Must read COMPOSITION before ITERATE"
exit op_loop
else if (.not.interaction_flag) then
write(6,*) "Error: Must read INTERACTION before ITERATE"
exit op_loop
else if (.not.unit_cell_flag) then
write(6,*) "Error: Must read UNIT_CELL before ITERATE"
exit op_loop
else if (.not.discretize_flag) then
write(6,*) "Error: Must read DISCRETIZATION before ITERATE"
exit op_loop
else if (.not.basis_flag) then
write(6,*) "Error: Must read BASIS before ITERATE"
exit op_loop
end if
! Read omega file, if not input preceding RESCALE command
if (.not.omega_flag) then
call input(input_filename, 'input_filename') ! input file prefix
open(unit=field_unit,file=trim(input_filename),&
status='old',iostat=ierr)
if (ierr/=0) stop "Error while opening omega source file."
call input_field(omega,field_unit)
close(field_unit)
omega_flag = .TRUE.
end if
call input(output_prefix,'output_prefix') ! output file prefix
iterate_flag = .TRUE.
call cpu_time(basis_time)
basis_time = basis_time - start_time
call cpu_time(start_time)
! Read parameters for iteration
call input_iterate_param
! Allocate and initialize chain objects used in scf_mod
! Create fft_plan, which is saved in scf_mod as public variable
call density_startup(ngrid, extr_order, chain_step,&
update_chain=.false.)
if (itr_algo=='NR') then
! Allocate private arrays for Newton-Raphson iteration
call iterate_NR_startup(N_star)
write(6,FMT = "( / '************************************' / )" )
! Main Newton-Raphson iteration loop
call iterate_NR( &
N_star, &! # of basis functions
omega, &! chemical potential field (IN/OUT)
itr, &! actual number of interations
converge, &! = .TRUE. if converged
error, &! final error = max(residuals)
rho, &! monomer density field
f_Helmholtz, &! Helmholtz free energy per monomer/kT
pressure, &! pressure * monomer volume / kT
stress &! d(free energy)/d(cell parameters)
)
else if(itr_algo=='AM') then
! Allocate private arrays for Anderson-Mixing iteration
call iterate_AM_startup(N_star)
call iterate_AM( &
N_star, &! # of basis functions
omega, &! chemical potential field (IN/OUT)
itr, &! actual number of interations
converge, &! = .TRUE. if converged
error, &! final error = max(residuals)
rho, &! monomer density field
f_Helmholtz, &! Helmholtz free energy per monomer/kT
pressure, &! pressure * monomer volume / kT
stress &! d(free energy)/d(cell parameters)
)
endif
! Defer output to beginning of next operation
! If next operation is SWEEP, '0.' will be added to output_prefix
output_flag = .TRUE.
call cpu_time(scf_time)
scf_time = scf_time - start_time
case ('SWEEP')
! Iterate to convergences for a sequence of sets of parameters.
! See sweep_mod.
! Must be preceded by successful ITERATE for first solution.
if (.not.iterate_flag) then
write(6,*) &
"Error: Must call ITERATE (1st iteration) before SWEEP"
exit op_loop
else if (.not.converge) then
write(6,*) "Error: 1st iteration failed to converge"
exit op_loop
end if
sweep_flag = .TRUE.
! Read parameters needed by sweep
call input(s_max,'s_max') ! max(contour variable s)
call input_increments(5,'N',domain) ! see sweep_mod
! Initialize contour variable s = 0.0 -> s_max
s = 0.0_long
! Initialize history arrays
call history_setup
call update_history(s,omega,cell_param,domain)
! Loop over sweep through parameters
i = 0
step_unit = 1.0_long
sweep_loop : do
if (i == 0) then
step = 0.1*step_unit
else if (i == 1) then
if (converge) then
step = 0.9*step_unit
else
step = step_unit - s
end if
else if (i > 1) then
step = step_unit
end if
s = s + step
call increment_parameters(step, domain, cell_param)
write(6, FMT = "('************************************'/ )" )
write(6, FMT = "('s =',f10.4)" ) s
call cpu_time(start_time)
! 1st order continuation of omega and cell_param
call continuation(step, domain, omega, cell_param)
! Reconstruct unit cell and all values of |k|^2
call make_unit_cell
call make_ksq(G_basis)
! Rebuild chains
call density_startup(ngrid, extr_order, chain_step, &
update_chain=.TRUE.)
! Main iteration routine
if (itr_algo=='NR')then
call iterate_NR( &
N_star, &! # of basis functions
omega, &! chemical potential field (IN/OUT)
itr, &! actual number of interations
converge, &! = .TRUE. if converged
error, &! final error = max(residuals)
rho, &! monomer density field
f_Helmholtz, &! Helmholtz free energy per monomer/kT
pressure, &! pressure * monomer volume/kT
stress &! d(free energy)/d(cell parameters)
)
else if (itr_algo=='AM')then
call iterate_AM( &
N_star, &! # of basis functions
omega, &! chemical potential field (IN/OUT)
itr, &! actual number of interations
converge, &! = .TRUE. if converged
error, &! final error = max(residuals)
rho, &! monomer density field
f_Helmholtz, &! Helmholtz free energy per monomer/kT
pressure, &! pressure * monomer volume / kT
stress &! d(free energy)/d(cell parameters)
)
end if
if (converge) then
i = i + 1
call update_history(s, omega, cell_param, domain)
call cpu_time(scf_time)
scf_time = scf_time - start_time
! Output out, rho, and omega files if s is an integer
if ( abs(s-float(nint(s))) < 0.001_long ) then
j = nint(s)
! Appending 'j.' to output_prefix in file names
call output_summary( &
trim(output_prefix)//trim(int_string(j))//'.' )
call output_fields( &
trim(output_prefix)//trim(int_string(j))//'.' )
end if
write(6,*)
else if (step_unit > ( (1.0/16.0) + 0.001 ) ) then
! Backtrack to previous state point
s = s - step
call increment_parameters(-step, domain, cell_param)
! Halve step size
step_unit = 0.5*step_unit
write(6, FMT="( / 'Backtrack and halve step_unit' / )" )
cycle
else ! If not.converge and step_unit <= 1/16, stop.
write(6,*)
write(6,*) 'Failed to converge - stop program'
stop
end if
if ( (s + 0.001) >= s_max ) exit sweep_loop
end do sweep_loop
! end loop over sweep through parameters
case ('RESPONSE')
! Calculate and diagonalize linear SCF response functions for
! a periodic microstructure. See response_mod.
! Check preconditions (Must iterate to convergence first)
if (.not. iterate_flag) then
write(6,*) "Error: Must ITERATE before calculating RESPONSE"
exit op_loop
end if
iterate_flag = .FALSE.
call response_startup(Ngrid, chain_step, order=1)
call response_sweep(Ngrid, output_prefix)
call release_basis()
call make_basis(R_basis,G_basis,group_name,Ngrid,grid_flag=.TRUE.)
case ('OUTPUT_GROUP')
! Check preconditions (Needs group created in BASIS block)
if (.not. basis_flag) then
write(6,*) "Error: Must read BASIS before OUTPUT_GROUP"
exit op_loop
end if
iterate_flag = .FALSE.
omega_flag = .FALSE.
open(unit=field_unit, &
file=trim(output_prefix)//'group',status='replace')
call output_group(group,field_unit)
close(field_unit)
case ('OUTPUT_WAVES')
! Check preconditions (Needs BASIS)
if (.not. basis_flag) then
write(6,*) "Error: Must read BASIS before OUTPUT_WAVES"
exit op_loop
end if
open(unit=field_unit,file=trim(output_prefix)//'waves',&
status='replace')
call output_waves(field_unit, group_name)
close(field_unit)
case ('FIELD_TO_RGRID')
! Read representation of field as list of coefficients of
! symmetrized basis functions, output file containing field
! values at grid points (rgrid).
! Check preconditions for FIELD_TO_RGRID
if ( .not. unit_cell_flag ) then
write(6,*) "Error: Must read UNIT_CELL before FIELD_TO_GRID"
exit op_loop
else if (.not.discretize_flag) then
write(6,*) &
"Error: Must read DISCRETIZATION before FIELD_TO_GRID"
exit op_loop
else if (.not.basis_flag) then
write(6,*) "Error: Must read BASIS before FIELD_TO_GRID"
exit op_loop
end if
iterate_flag = .FALSE.
omega_flag = .FALSE.
! Read input and output file names from input script
call input(input_filename,'input_filename')
call input(output_filename,'output_filename')
! Read field (coefficients of basis functions) from input_filename
open(unit=field_unit,file=trim(input_filename),status='old')
call input_field(rho,field_unit)
close(field_unit)
! Write values of field on a grid to output_filename
open(unit=field_unit,file=trim(output_filename),status='replace')
call output_field_grid(rho,field_unit,group_name,ngrid)
close(field_unit)
case ('KGRID_TO_RGRID')
! Read representation of field as list of coefficients of
! on FFT grid (kgrid) and output file containing field
! values at grid points (rgrid).
! Check preconditions for KGRID_TO_RGRID
if ( .not. unit_cell_flag ) then
write(6,*) "Error: Must read UNIT_CELL before KGRID_TO_RGRID"
exit op_loop
else if (.not.discretize_flag) then
write(6,*) &
"Error: Must read DISCRETIZATION before KGRID_TO_RGRID"
exit op_loop
else if (.not.basis_flag) then
write(6,*) "Error: Must read BASIS before KGRID_TO_RGRID"
exit op_loop
end if
iterate_flag = .FALSE.
omega_flag = .FALSE.
! Read input and output file names from input script
call input(input_filename,'input_filename')
call input(output_filename,'output_filename')
if (.not.allocated(k_grid)) then
allocate(k_grid(0:ngrid(1)/2, 0:ngrid(2)-1, 0:ngrid(3)-1, N_monomer))
end if
! Open input_filename, read header, check grid dimensions
open(unit=field_unit,file=trim(input_filename),status='old')
! Skip first 13 lines
do i=1,14
read(field_unit,*)
end do
read (field_unit,*) grid_size
if (grid_size(1) /= ngrid(1)) then
write(6,*) "Error: Inconsistent grid in input kgrid file"
exit op_loop
end if
if (dim > 1) then
if (grid_size(2) /= ngrid(2)) then
write(6,*) "Error: Inconsistent grid in input kgrid file"
exit op_loop
end if
if (dim > 2) then
if (grid_size(3) /= ngrid(3)) then
write(6,*) "Error: Inconsistent grid in input kgrid file"
exit op_loop
end if
end if
end if
! Read elements of k_grid (Fourier coefficients) from input file
k_grid = 0.0
do i1 = 0, ngrid(1)/2
do i2 = 0, ngrid(2) - 1
do i3 = 0, ngrid(3) - 1
read(field_unit,*) k_grid(i1,i2,i3,:)
end do
end do
end do
close(field_unit)
! Transform to symmetry-adapated Fourier expansion of rho
call create_fft_plan(ngrid, plan)
do alpha=1, N_monomer
call kgrid_to_basis(k_grid(:,:,:,alpha), rho(alpha,:))
end do
! Write rho field in coordinate grid format
open(unit=field_unit,file=trim(output_filename),status='replace')
call output_field_grid(rho, field_unit, group_name, ngrid)
close(field_unit)
deallocate(k_grid)
case ('RGRID_TO_FIELD')
! Read coordinate-grid representation of a field
! Write representation in symmetry-adapated Fourier expansion
! Check preconditions for RGRID_TO_FIELD
if ( .not. unit_cell_flag ) then
write(6,*) "Error: Must read UNIT_CELL before RGRID_TO_FIELD"
exit op_loop
else if (.not.discretize_flag) then
write(6,*) &
"Error: Must read DISCRETIZATION before RGRID_TO_FIELD"
exit op_loop
else if (.not.basis_flag) then
write(6,*) "Error: Must read BASIS before RGRID_TO_FIELD"
exit op_loop
end if
iterate_flag = .FALSE.
omega_flag = .FALSE.
! Read input and output file names from input script
call input(input_filename,'input_filename')
call input(output_filename,'output_filename')
! Allocate required memory
if (.not.allocated(r_grid)) then
allocate(r_grid(0:ngrid(1)-1,0:ngrid(2)-1,0:ngrid(3)-1,N_monomer))
end if
if (.not.allocated(k_grid)) then
allocate(k_grid(0:ngrid(1)/2,0:ngrid(2)-1,0:ngrid(3)-1,N_monomer))
end if
! Open input file and read header, including grid dimensions
open(unit=field_unit,file=trim(input_filename),status='old')
! Skip first 13 lines
do i=1,14
read(field_unit, *)
end do
read(field_unit,*) grid_size
if (grid_size(1) /= ngrid(1)) then
write(6,*) "Error: inconsistent grid size in kgrid input file"
exit op_loop
endif
if (dim > 1) then
if (grid_size(2) /= ngrid(2)) then
write(6,*) "Error: Inconsistent grid size in rho_kgrid file"
exit op_loop
end if
if (dim > 2) then
if (grid_size(3) /= ngrid(3)) then
write(6,*) "Error: Inconsistent grid size in kgrid file"
exit op_loop
end if
end if
end if
! Read field values at grid points from input file
r_grid=0.0
do i3 = 0, ngrid(3) - 1
do i2 = 0, ngrid(2) - 1
do i1 = 0, ngrid(1) - 1
read(field_unit,*) r_grid(i1,i2,i3,:)
end do
end do
end do
close(field_unit)
! Transform to symmetry-adapated rho field
call create_fft_plan(ngrid, plan)
rnodes=dble( plan%n(1) * plan%n(2) * plan%n(3) )
do alpha = 1, N_monomer
call fft(plan,r_grid(:,:,:,alpha), k_grid(:,:,:,alpha))
k_grid(:,:,:,alpha) = k_grid(:,:,:,alpha)/rnodes
call kgrid_to_basis(k_grid(:,:,:,alpha), rho(alpha,:))
end do
open(unit=field_unit, file=trim(output_filename), status='replace')
call output_field(rho, field_unit, group_name)
close(field_unit)
deallocate(r_grid)
deallocate(k_grid)
case('RHO_TO_OMEGA')
! (1) Read a rho field from file, in symmetry-adapted format.
! (2) Generate approximate omega field from rho field by assuming
! that Lagrange multiplier (pressure) field is zero.
! (3) Output resulting omega file to file in symmetry-adapated format
! Check preconditions for RGRID_TO_FIELD
if ( .not. unit_cell_flag ) then
write(6,*) "Error: Must read UNIT_CELL before RHO_TO_OMEGA"
exit op_loop
else if (.not.discretize_flag) then
write(6,*) &
"Error: Must read DISCRETIZATION before RHO_TO_OMEGA"
exit op_loop
else if (.not.basis_flag) then
write(6,*) "Error: Must read BASIS before RHO_TO_OMEGA"
exit op_loop
end if
iterate_flag = .FALSE.
omega_flag = .FALSE.
! Read filenames from parameter file
call input(input_filename, 'input_filename')
call input(output_filename, 'output_filename')
! Read input rho field
open(unit=field_unit,file=trim(input_filename),status='old')
call input_field(rho,field_unit)
close(field_unit)
! Compute approximate omega field
allocate(omega_basis(N_monomer, N_star))
do alpha=1, N_monomer
do i=1, N_star
omega_basis(alpha,i) = sum(chi(:,alpha)*rho(:,i))
end do
end do
! Output omega field
open(unit=field_unit,file=trim(output_filename),status='replace')
call output_field(omega_basis,field_unit,group_name)
close(field_unit)
deallocate(omega_basis)
case ('FINISH')
! Stop execution
exit op_loop
case default
write(6,*) 'Error: Invalid op_string'
exit op_loop
end select
end do op_loop
contains ! internal subroutines of program pscf_pd
subroutine allocate_scf_arrays
!-----------------------------------------------------------------
! Allocate arrays omega, rho, stress
!-----------------------------------------------------------------
integer :: i ! error for file opening
allocate(omega(N_monomer,N_star), stat = i)
if (i.ne.0 ) stop 'Error allocating omega'
allocate(rho(N_monomer,N_star), stat = i)
if (i.ne.0) stop 'Error allocating rho'
allocate(stress(N_cell_param), stat = i )
if (i.ne.0) stop 'Error allocating stress'
end subroutine allocate_scf_arrays
!=================================================================
subroutine output_summary(prefix)
!-----------------------------------------------------------------
! Writes the output summary to the file output_prefix//suffix
!-----------------------------------------------------------------
use chemistry_mod, only : output_chemistry, N_chain, N_solvent, &
ensemble,phi_chain,phi_solvent,mu_chain, &
mu_solvent, interaction_type, chi
use scf_mod, only : free_energy_FH
character(*) :: prefix
real(long) :: f_homo ! FH free energy, kT / monomer
open(file=trim(prefix)//'out',unit=out_unit,status='replace')
call set_io_units(o=out_unit)
! Set version number for *.out file
version%major = 1
version%minor = 0
call output_version(version, out_unit)
! Output in format of input driver file
if (monomer_flag) then
write(out_unit,*)
call output('MONOMERS',f='N',j='L')
call output_monomers(out_unit,'F')
end if
if (chain_flag) then
write(out_unit,*)
call output('CHAINS',f='N',j='L')
call output_chains(out_unit,'F')
end if
if (solvent_flag) then
write(out_unit,*)
call output('SOLVENTS',f='N',j='L')
call output_solvents(out_unit,'F')
end if
if (composition_flag) then
write(out_unit,*)
call output('COMPOSITION',f='N',j='L')
call output_composition(out_unit,'F')
end if
if (interaction_flag) then
write(out_unit,*)
call output('INTERACTION',f='N',j='L')
call output_interaction(out_unit,'F')
end if
if (unit_cell_flag) then
write(out_unit,*)
call output('UNIT_CELL',f='N',j='L')
call output_unit_cell(out_unit,'F')
end if
if (discretize_flag) then
write(out_unit,*)
call output('DISCRETIZATION',f='N',j='L')
call output(ngrid,dim,'ngrid')
call output(chain_step,'chain_step')
!call output(extr_order, 'extr_order')
end if
if (basis_flag) then
write(out_unit,*)
call output('BASIS',f='N',j='L')
call output(trim(group_name),'group_name')
end if
if (iterate_flag) then
write(out_unit,*)
call output('ITERATE',f='N',j='L')
call output(trim(input_filename), 'input_filename')
call output(trim(output_prefix),'output_prefix')
call output_iterate_param
end if
write(out_unit,*)
call output('FINISH',f='N',j='L')
! End input script section, begin additional information
! Thermodynamics
write(out_unit,*)
call output('THERMO',f='N',j='L')
call output(f_Helmholtz,'f_Helmholtz')
f_homo = free_energy_FH(phi_chain,phi_solvent)
call output(f_homo, 'f_homo')
call output(pressure, 'pressure')
select case(ensemble)
case(0) ! phi was output by output chemistry, so
if( N_chain > 0 ) then
call output(mu_chain,N_chain,'mu_chain',s='C')
end if
if( N_solvent > 0 ) then
call output(mu_solvent,N_solvent,'mu_solvent',s='C')
end if
case(1) ! mu was output chemistry, so
if( N_chain > 0) then
call output(phi_chain,N_chain,'phi_chain',s='C')
end if
if ( N_solvent > 0) then
call output(phi_solvent,N_solvent,'phi_solvent',s='C')
end if
end select
call output(stress,N_cell_param,'stress')
! Output chi if chi is input as chi = chi_A/T + B
if (interaction_type =='chi_T') then
call output(chi,N_monomer,N_monomer,'chi',s='L')
end if
! Timing and resource statistics
write(out_unit,*)
call output('STATISTICS',f='N',j='L')
call output(N_star,'N_star')
call output(error,'Final Error')
call output(itr,'Iterations')
call output(basis_time,'Basis Time')
call output(scf_time,'SCF Time')
close(out_unit) ! close output_prefix//'out' file
call set_io_units(o=6) ! reset default echo unit to stdout
end subroutine output_summary
!============================================================
subroutine output_fields(prefix)
!-----------------------------------------------------------------
! Writes the output summary to the file output_prefix//suffix
!-----------------------------------------------------------------
character(*) :: prefix
open(file=trim(prefix)//'omega', unit=field_unit,status='replace')
call output_field(omega,field_unit,group_name)
close(field_unit)
open(file=trim(prefix)//'rho',unit=field_unit,status='replace')
call output_field(rho,field_unit,group_name)
close(field_unit)
end subroutine output_fields
!============================================================
end program pscf_pd