[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
!---------------------------------------------------------------------