[table of contents] [master index] [comments] [modules] [programs] [variables] [types] [procedures]
MODULE
iterate_mod
PURPOSE
Iteration of SCF equations
ALGORITHM
Newton/Broyden algorithm that uses an approximate calculation for the initial Jacobian and Broyden updates to the inverse Jacobian. The approximate initial Jacobian is calculatd from a corresponding approximate linear response matrix, which is calculated in response_pd_mod. See the docblock for N_cut for a discussion of approximate response function. Subroutine iterate_NR can either: i) Iterate omega in a fixed unit cell (domain = F) ii) Iterate both omega and unit_cell_param (domain = T) depending upon value of logical module variable domain
SOURCE
module iterate_mod use const_mod use io_mod use response_pd_mod implicit none private ! public procedures public :: input_iterate_param ! input iteration parameters public :: output_iterate_param ! output iteration variables public :: iterate_NR_startup ! allocates arrays needed by iterate_NR public :: iterate_AM_startup ! allocates arrays needed by iterate_AM public :: iterate_NR ! main Newton-Raphson iteration routine public :: iterate_AM ! main Anderson-Mixing iteration routine ! public variables public :: max_itr ! maximum allowable number of iterations public :: error_max ! error tolerance: stop when error < error_max public :: domain ! If true, adjust unit_cell_dimensions public :: itr_algo ! Iteration algorithm public :: N_cut ! dimension of cutoff response and Jacobian for NR scheme public :: N_hist ! Number of previous histories to use in Anderson-Mixing Scheme
VARIABLE
integer max_itr - maximum allowed number of iterations
VARIABLE
real(long) error_max - maximum error in converged solution - Stop when error < error_max
VARIABLE
logical domain - If .true., iterate variable unit_cell shape If .false., iterate omega in fixed unit cell
VARIABLE
character(10) iter_algo - Iteration algorithm to be used For now, must equal 'NR' (Newton-Raphson)
PURPOSE
This variable will eventually allow the user to select from more than one iteration algorithm. Thus far, however, the only algorithm that has been implemented is a quasi-Newton / Broyden algorithm, for which the value of iter_algo should be 'NR'. Implementation of a more robust relaxation algorithm is planned.
VARIABLE
integer N_cut - dimension of cutoff response and Jacobian
PURPOSE
The approximate linear response matrix used to construct the initial Jacobian is obtained by numerically calculating the response to the first N_cut basis functions to obtain an N_cut x N_cut submatrix within the subspace spanned by these functions, and using the analytically calculated response of a homogenous gas to approximate the response to the remaining basis functions.
SUBROUTINE
input_iterate_param
PURPOSE
Read parameters max_itr, error_max, domain, itr_algo, N_cut
SOURCE
subroutine input_iterate_param
SUBROUTINE
output_iterate_param
PURPOSE
Read parameters max_itr, error_max, domain, N_cut
SOURCE
subroutine output_iterate_param
SUBROUTINE
iterate_NR_startup(N)
PURPOSE
1) Set integer module variable N_residual = # of residuals 2) Allocate arrays needed by iterate_NR
ARGUMENTS
integer N = # of basis functions
SOURCE
subroutine iterate_NR_startup(N) use unit_cell_mod, only : N_cell_param use chemistry_mod, only : N_monomer, ensemble integer, intent(IN) :: N ! # of basis functions
SUBROUTINE
subroutine iterate_NR
PURPOSE
Newton-Raphson iteration of the SCFT, using numerically calculated Jacobian
COMMENT
1) Algorithm assumes that density_startup, iterate_NR_startup, and make_unit_cell have been called prior to iterate_NR. 2) Routine will simultaneously iterates omega field and unit cell parameters if domain = .true., or iterate the omega field for a fixed unit cell if domain = .false.
SOURCE
subroutine iterate_NR( & N, &! # 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) ) !-------------------------------------------------------------------- use io_mod, only : output use unit_cell_mod, only : N_cell_param, cell_param, & make_unit_cell, G_basis, dGG_basis use basis_mod, only : make_dGsq use chemistry_mod, only : N_monomer, ensemble, mu_chain, mu_solvent, & phi_solvent, phi_chain, N_chain, N_solvent use scf_mod, only : density, scf_stress, set_omega_uniform, & mu_phi_chain, mu_phi_solvent, free_energy use grid_mod, only : make_ksq ! Arguments integer, intent(IN) :: N real(long), intent(INOUT) :: omega(:,:) ! omega(N_monomer,N) integer, intent(OUT) :: itr logical, intent(OUT) :: converge real(long), intent(OUT) :: error real(long), intent(OUT) :: rho(:,:) ! rho(N_monomer,N) real(long), intent(OUT) :: f_Helmholtz real(long), intent(OUT) :: pressure real(long), intent(OUT) :: stress(:) ! stress(N_cell_param)
SUBROUTINE
make_residual( N, domain, rho, stress, omega, residual )
PURPOSE
Calculate array residual(1:N_residual) Uses input values of rho, omega, and stress arrays
ARGUMENTS
integer N = # of basis functions logical domain = T -> domain iteration, F -> fixed cell real(long) rho(:,:) = monomer density field real(long) stress(:) = stress components real(long) omega(:,:) = chemical potential field real(long) residual(:) = residual array, dimension(1:N_residual)
SOURCE
subroutine make_residual( N, domain, rho, stress, omega, residual ) use unit_cell_mod, only: N_cell_param use chemistry_mod, only: ensemble, N_monomer, chi implicit none integer, intent(IN) :: N logical, intent(IN) :: domain real(long), intent(IN) :: rho(:,:) real(long), intent(IN) :: stress(:) real(long), intent(IN) :: omega(:,:) real(long), intent(OUT):: residual(:)
SUBROUTINE
Jacobian_response(N, N_cut, domain, omega, Jacobian)
PURPOSE
Construct Jacobian from response functions
SOURCE
subroutine Jacobian_response(N, N_cut, domain, omega, Jacobian) use chemistry_mod, only : ensemble, N_monomer, chi use unit_cell_mod, only : N_cell_param implicit none integer, intent(IN) :: N integer, intent(IN) :: N_cut logical, intent(IN) :: domain real(long), intent(IN) :: omega(:,:) real(long), intent(OUT) :: Jacobian(:,:)
SUBROUTINE
iterate_AM_startup(N)
PURPOSE
1) Allocate arrays needed by iterate_AM
ARGUMENTS
integer N = # of basis functions
SOURCE
subroutine iterate_AM_startup(N) use unit_cell_mod, only : N_cell_param use chemistry_mod, only : N_monomer, ensemble integer, intent(IN) :: N ! # of basis functions integer :: info ! message variable for dgetri !! real(long), allocatable :: Umn(:,:) ! Local Umn variable to calculate the ! work space if (allocated(ipvt)) deallocate(ipvt) if (allocated(work)) deallocate(work) if (allocated(dev)) deallocate(dev) if (allocated(omega_hist)) deallocate(omega_hist) if (allocated(Umn)) deallocate(Umn) allocate(ipvt(N_hist)) allocate(work(N_hist)) allocate(dev(N_hist+1,N-1,N_monomer)) allocate(omega_hist(N_hist+1,N-1,N_monomer)) allocate(Umn(N_hist,N_hist)) ! estimate the opitimal workspace dimension lwork = -1 call dgetri(N_hist,Umn,N_hist,ipvt,work,lwork,info) lwork = work(1) if (allocated(work)) deallocate(work) allocate(work(lwork)) ! Umn changes size after every iteration so will be alloated in main (next) ! subroutine, however, here we allocated it to calculate the maximum work ! space needed for the AM scheme. if (allocated(Umn)) deallocate(Umn) end subroutine iterate_AM_startup !---------------------------------------------------------------------