From 129f6e4625d150493e127dee6a110e0241aa4e3c Mon Sep 17 00:00:00 2001 From: David Morse Date: Tue, 14 Nov 2017 21:37:27 -0600 Subject: [PATCH 1/9] Editing README.md [ci skip] --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 66193ed..503b156 100644 --- a/README.md +++ b/README.md @@ -56,7 +56,7 @@ If you use PSCF in published work, we request that you cite the paper: "Broadly Accessible Self-Consistent Field Theory for Block Copolymer Materials Discovery", A. Arora, J. Qin, D.C. Morse, K.T. Delaney, G.H. Fredrickson, F.S. Bates and K.D. Dorfman, -Macromolecules 49, 4675 (2016) +*Macromolecules* **49**, 4675 (2016) available electronically at http://pubs.acs.org/doi/10.1021/acs.macromol.6b00107 From 8059761ce3b111cdcba6be733fe3e192dfafe007 Mon Sep 17 00:00:00 2001 From: David Morse Date: Tue, 14 Nov 2017 21:41:08 -0600 Subject: [PATCH 2/9] Editing README.md [ci skip] --- README.md | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/README.md b/README.md index 503b156..14b0d44 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # PSCF - Polymer Self-Consistent Field Theory -Copyright (2002-2016) Regents of the University of Minnesota +Copyright (2002-2017) Regents of the University of Minnesota PSCF is a Fortran 90 program for numerically solving the polymer self-consistent field theory (SCFT) for spatially periodic structures @@ -12,15 +12,6 @@ it under the terms of the GNU General Public License as published by the Free Software Foundation. A copy of this license is included in the LICENSE file in the top-level PSCF directory. -# Contributors - -- David Morse -- Chris Tyler -- Jian Qin -- Amit Ranjan -- Raghuram Thiagarajan -- Akash Arora - # Dependencies PSCF depends upon the FFTW fast Fourier transform library and the @@ -91,3 +82,12 @@ An annotated list of source files is given in the file src/SRC_FILES. Before modifying any fortran files, see the note at the end of that file regarding the use of preprocessor to generate some files. +# Contributors + +- David Morse +- Chris Tyler +- Jian Qin +- Amit Ranjan +- Raghuram Thiagarajan +- Akash Arora + From 71538554bb32b61996832a020a287853edd324b0 Mon Sep 17 00:00:00 2001 From: David Morse Date: Tue, 14 Nov 2017 22:00:56 -0600 Subject: [PATCH 3/9] Editing param.rst manual page --- doc/user-man/param.rst | 66 ++++++++++++++++++++++-------------------- 1 file changed, 35 insertions(+), 31 deletions(-) diff --git a/doc/user-man/param.rst b/doc/user-man/param.rst index 39ac98d..21385d7 100644 --- a/doc/user-man/param.rst +++ b/doc/user-man/param.rst @@ -6,15 +6,17 @@ Parameter File ************** The main program reads an parameter file containing the parameters and -instructions for a calculation. This file is divided into sections, -each of which contains a different type of information. Each section -is preceded by a blank line and starts with a line containing a -section title in all capital letters (i.e., 'CHEMISTRY', 'UNIT_CELL', -etc.) Each block may contain values of a sequence of variables. The -name of each variable appears on a line by itself, followed by the -value or (for arrays) values on one more subsequent lines. The -order in which variables must appear within a section is fixed. The -program stops when it encounters the block title 'FINISH'. +instructions for a calculation. This file is divided into sections, +each of which contains a different type of information. + +Each section is preceded by a blank line and starts with a line containing +a section title in all capital letters (i.e., 'CHEMISTRY', 'UNIT_CELL', +etc.) Each section may contain values of a sequence of variables. The +name of each variable appears on a line by itself, followed by the value +or (for arrays) values on one more subsequent lines. The variables within +each section must appear in a predetermined (i.e., hard-coded) order. + +The program stops when it encounters the section title 'FINISH'. .. _example-sec: @@ -28,16 +30,15 @@ the blocks. The first line of the file identifies the version of the file format (in this case, version 1.0). The remainder of the file is divided into -sections, each of which begins with a line containing a capitalized label, -such as MONOMERS, CHAINS, ec. Each section begins with a capitalized -section identtifier (MONOMER, CHAINS, etc.) on a line by itself. A single -line must is left blank between sections. The sections are processed in -the order in which the appear in the parameter file. The first few sections -in this example simply provide values for physical and computational -parameters. An ITERATE section instructs the program to actually perform -a SCF calculation, by iteratively solving the SCF equations. A SWEEP -section performs a sequence of such calculations along a line in parameter -space. Execution of the program stops when a FINISH line is encountered. +sections. Each section begins with a capitalized section identtifier +(MONOMER, CHAINS, etc.) on a line by itself. A single blank line appears +between sections. The sections are processed in the order in which they +appear in the parameter file. The first few sections in this example +simply provide values for physical and computational parameters. An +ITERATE section instructs the program to actually perform a single SCF +calculation, by iteratively solving the SCF equations. A SWEEP section +performs a sequence of such calculations along a line in parameter space. +Execution of the program stops when a FINISH line is encountered. :: @@ -212,20 +213,23 @@ or more of the statistical segment lengths is set to unity. SCFT also leaves the user some freedom to redefine what he or she means by a "monomer", which need not correspond to a chemical repeat unit. The choice of values of the parameters block_length, solvent_size, kuhn, and - chi to represent a particular experimental system all depend on the choice -of a value for a monomer reference volume, which in turn defines an effective -monomer repeat unit. A monomer of a polymeric species defined to be a length -or molar mass of chain that occupies one monomer reference volume in the melt. -Each element of the variable block_length represents the number of "monomers" -in a block of a block copolymer, defined to be the ratio of the block volume -to the monomer reference volume. Similarly, the variable solvent_size is given -by ratio of the solvent volume to the reference volume. The values of the chi -parameters are proportional to the reference volume, while kuhn lengths are -proportional to the square root of the reference volume. +chi to represent a particular experimental system all depend on an implicit +choice of a value for a monomer reference volume, which defines the mononmer +repeat unit that is being used for bookkeeping purposes. One "monomer" of a +polymeric species is defined to correspond to length or molar mass of chain +that occupies a volume in the melt equal to one reference volume, which may +or may not correspond to a chemical repeat unit. Each element of the variable +block_length represents the number of "monomers" in a block of a block +copolymer, which is given by the ratio of the block volume to the monomer +reference volume. Similarly, the variable solvent_size is given by ratio +of the solvent volume to the reference volume. Note that PSCF does not require the user to input a value for the monomer -reference volume - the choice is implicit in the values given for other -quantities. +reference volume - the choice of reference volume is implicit in the values +given for other quantities. Changes in ones choice of reference volume lead +to corresponding changes in the values for the chi parameters, which are +proportional to the reference volume, and in the kuhn lengths, which are +proportional to the square root of the reference volume. **Character Strings** From c57822380f709139de1d2a9643636edea9602af5 Mon Sep 17 00:00:00 2001 From: David Morse Date: Tue, 14 Nov 2017 22:06:55 -0600 Subject: [PATCH 4/9] Editing field.rst [ci skip] --- doc/user-man/field.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/user-man/field.rst b/doc/user-man/field.rst index d254272..f84fe87 100644 --- a/doc/user-man/field.rst +++ b/doc/user-man/field.rst @@ -52,11 +52,11 @@ the field :math:`\omega_{\alpha}` as an expansion of the form .. math:: \omega_{\alpha}(\textbf{r}) = - \sum_{i=1}^{\textrm{N\_star}} \omega_{i\alpha} f_{i}(\textbf{r}) + \sum_{i=1}^{N_{star}} \omega_{i\alpha} f_{i}(\textbf{r}) in which each function :math:`f_{i}(\textbf{r})` is a real basis function, :math:`\omega_{i\alpha}` is an associated real coefficient, -and :math:`\textrm{N\_star}` is the number of basis functions used to +and :math:`N_{star}` is the number of basis functions used to approximate the field. In a symmetry-adapted Fourier expansion of a field with a specified space group symmetry, each basis function :math:`f_{i}(\textbf{r})` is a real function that is invariant under From f50a19641befbced4f99a552972e9e63d68a58ac Mon Sep 17 00:00:00 2001 From: David Morse Date: Tue, 14 Nov 2017 22:11:03 -0600 Subject: [PATCH 5/9] Editing index.rst [ci skip] --- doc/user-man/index.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/user-man/index.rst b/doc/user-man/index.rst index d73ee39..de74ffe 100644 --- a/doc/user-man/index.rst +++ b/doc/user-man/index.rst @@ -3,8 +3,8 @@ You can adapt this file completely to your liking, but it should at least contain the root `toctree` directive. -Welcome to PSCF's documentation! -================================ +PSCF User Manual +================ Contents: From 5d4ea149a2ebf30ed2ec10df049b8bb4d948bae4 Mon Sep 17 00:00:00 2001 From: "J. Haataja" Date: Wed, 13 Dec 2017 17:08:57 +0200 Subject: [PATCH 6/9] Fix typo in linux installation instructions + additional recommended package. --- doc/user-man/binary.rst | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/doc/user-man/binary.rst b/doc/user-man/binary.rst index 91b8faa..08e2875 100644 --- a/doc/user-man/binary.rst +++ b/doc/user-man/binary.rst @@ -70,11 +70,13 @@ files. To install from binary on Ubuntu: - libfftw3-3 - liblapack3 + - libatlas-base-dev To install these packages using apt-get, enter:: - sudo apt-get libfftw3-3 - sudo apt-get liblapack3 + sudo apt-get install libfftw3-3 + sudo apt-get install liblapack3 + sudo apt-get install libatlas-base-dev * Download the .deb package from the PSCF home page. This is a file with a name of the form pscf--Linux.deb, where From 48e9b0efb3a9e274797d032772fae93a69ae5761 Mon Sep 17 00:00:00 2001 From: akasharora123 Date: Sat, 2 Jun 2018 10:39:24 -0500 Subject: [PATCH 7/9] Anderson-Mixing variable cell implemented. I/O module modified to display 40 char long comments during iterations. --- src/io/io_mod.f | 4 +- src/iterate/iterate_mod.fp.f | 648 +++++++++++++++++------------------ 2 files changed, 314 insertions(+), 338 deletions(-) diff --git a/src/io/io_mod.f b/src/io/io_mod.f index 834316f..b0f1424 100644 --- a/src/io/io_mod.f +++ b/src/io/io_mod.f @@ -356,7 +356,7 @@ module io_mod character(1) :: default_mat_sym = 'N' ! normal (full matrix) ! Integer field widths for output of comments and data: - integer :: com_width = 20 ! comment field width + integer :: com_width = 40 ! comment field width integer :: data_width = 20 ! data field width (total) integer :: frac_width = 12 ! # digits after decimal @@ -364,7 +364,7 @@ module io_mod character(2) :: fmt_ef = 'ES' ! must = F, E, or ES ! Output format strings for scalar variables: - character(3) :: fmt_c = 'A20' ! comment format + character(3) :: fmt_c = 'A40' ! comment format character(3) :: fmt_i = 'I20' ! integer format character(7) :: fmt_r = 'ES20.10' ! real format character(3) :: fmt_l = 'L20' ! logical format diff --git a/src/iterate/iterate_mod.fp.f b/src/iterate/iterate_mod.fp.f index f35e2ac..de2307a 100644 --- a/src/iterate/iterate_mod.fp.f +++ b/src/iterate/iterate_mod.fp.f @@ -14,15 +14,21 @@ ! 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. +! ALGORITHMS +! 1. 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. +! +! 2. Anderson Mixing (J. Chem. Phys. 146, 244902 (2017)) +! A mixing algorithm that uses results of a few previous iterations +! to produce the guess for next iteration. Convergence depends on +! two parameters: lamdda (mixing parameter), and N_hist (Number of +! iterations). ! -! Subroutine iterate_NR can either: +! Subroutine iterate_NR/iterate_AM 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 @@ -31,6 +37,7 @@ ! Amit Ranjan - multi-component blend version (2003) ! David Morse - modifications (1/2004) ! Jian Qin - Broyden and backtracking (2006-2007) +! Akash Arora - Anderson Mixing ! SOURCE !----------------------------------------------------------------------- module iterate_mod @@ -122,7 +129,9 @@ module iterate_mod ! Varaibles for Anderson Mixing scheme real(long),allocatable :: Umn(:,:) real(long),allocatable :: dev(:,:,:) - real(long),allocatable :: omega_hist(:,:,:) + real(long),allocatable :: dev_cp(:,:) + real(long),allocatable :: om_hist(:,:,:) + real(long),allocatable :: cp_hist(:,:) contains @@ -387,7 +396,9 @@ subroutine iterate_NR( & write(6,"('Iteration ',i3)") itr ! Output to log file - call output(error, 'error =',o=6,f='L') +!! call output(error, 'error =',o=6,f='L') + call output(error,'max (SCF + 100xStress) residuals =',o=6,f='L') + call output(error1,'max SCF residuals =',o=6,f='L') if (domain) then call output(stress,N_cell_param, 'stress =',o=6,f='L') call output(cell_param,N_cell_param,'cell_param =',o=6,f='L') @@ -981,13 +992,17 @@ subroutine iterate_AM_startup(N) if (allocated(ipvt)) deallocate(ipvt) if (allocated(work)) deallocate(work) if (allocated(dev)) deallocate(dev) - if (allocated(omega_hist)) deallocate(omega_hist) + if (allocated(dev_cp)) deallocate(dev_cp) + if (allocated(om_hist)) deallocate(om_hist) + if (allocated(cp_hist)) deallocate(cp_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(dev_cp(N_hist+1,N_cell_param)) + allocate(om_hist(N_hist+1,N-1,N_monomer)) + allocate(cp_hist(N_hist+1,N_cell_param)) allocate(Umn(N_hist,N_hist)) @@ -1075,7 +1090,7 @@ subroutine iterate_AM( & ! parameter real(long), parameter :: large = 100.0 - real(long), parameter :: stress_rescale=1.0d+2 + real(long), parameter :: stress_rescale=1.0 real(long), parameter :: STPMX=100.0_long ! local Variables @@ -1088,14 +1103,19 @@ subroutine iterate_AM( & ! Parameters for Anderson-Mixing real(long), allocatable :: Vm(:,:), Cn(:,:) - real(long) :: WW(N-1), DD(N-1) - real(long) :: elm(N-1) + real(long) :: Vm_2, Cn_2, Umn_2 + real(long) :: WW(N-1), DD(N-1), WW_cp(N_cell_param), DD_cp(N_cell_param) + real(long) :: elm(N-1), elm_cp(N_cell_param) real(long) :: dev_2, omega_2 - real(long) :: lam_AM + real(long) :: lam_AM, f_old + real(long) :: error1, error2, error_new, error_new_1, error_A, error_B, error_C integer :: mm,nn, itr_prev dev=0 - omega_hist=0 + om_hist=0 + dev_cp=0 + cp_hist=0 + M = N - 1 + ensemble @@ -1115,78 +1135,114 @@ subroutine iterate_AM( & if (ensemble == 0) call set_omega_uniform(omega) itr = 0 - - ! AM scheme needs atleast two previous iterations to mix and produce - ! output for the next. Here, we use first iteration as a guess supplied - ! and then create two (although need one) more iterations by linear - ! combinations of the previous iterations - write(6,*)'Creating History for Anderson-Mixing Scheme' - do i=1,3 + do alpha=1,N_monomer + om_hist(1,:,alpha) = omega(alpha,2:) + end do + if(domain) cp_hist(1,:) = cell_param(1:N_cell_param) + + iterative_loop : do + itr = itr + 1 - do alpha=1,N_monomer - omega(alpha,2:) = 0.9**(itr-1)*omega(alpha,2:) - end do + ! Calculating Lambda value need for mixing + if(itr 0) then + call mu_phi_solvent(mu_solvent, phi_solvent, q_solvent) + endif + call free_energy(N, rho, omega, phi_chain, mu_chain, & + phi_solvent, mu_solvent, f_Helmholtz, pressure) + !# ifdef DEVEL + call divide_energy(rho,omega,phi_chain,phi_solvent,Q,f_Helmholtz,f_component,overlap) + !# endif + + ! Calcuting Deviations do alpha=1,N_monomer if(N_monomer==2)then do k=1,N-1 - dev(itr,k,alpha) = sum(chi(:,alpha)*rho(:,k+1))&! + dev(itr_prev+1,k,alpha) = sum(chi(:,alpha)*rho(:,k+1))&! + sum(omega(:,k+1))/N_monomer - omega(alpha,k+1) end do - + elseif(N_monomer==3)then do k=1,N-1 - dev(itr,k,alpha) = sum(chi(:,alpha)*rho(:,k+1)) &! + dev(itr_prev+1,k,alpha) = sum(chi(:,alpha)*rho(:,k+1)) &! + ( sum(omega(:,k+1)) + chi(1,2)*rho(3,k+1) &! + chi(2,3)*rho(1,k+1) + chi(1,3)*rho(2,k+1) )/N_monomer &! - omega(alpha,k+1) end do endif - omega_hist(itr,:,alpha) = omega(alpha,2:) - - end do - end do - itr = itr - 1 ! i=1 in the previous loop is 0th iteration - - ! Calculating error in omega - dev_2 = 0.0 - omega_2 = 0.0 - do alpha = 1,N_monomer - dev_2 = dev_2 + norm_2(dev(itr,:,alpha))**2.0 - omega_2 = omega_2 + norm_2(omega(alpha,2:))**2.0 - end do - error = (dev_2/omega_2)**0.5 - ! Calculate thermodynamic properties - call free_energy(N,rho,omega,phi_chain,mu_chain,phi_solvent,mu_solvent,f_Helmholtz,pressure) - !# ifdef DEVEL - call divide_energy(rho,omega,phi_chain,phi_solvent,Q,f_Helmholtz,f_component,overlap) - !# endif - if (domain) stress = scf_stress(N, N_cell_param, dGsq) - + end do + + if(domain) dev_cp(itr_prev+1,:) = -stress - itr_stress = 0 - stress_loop : do + ! Calculating Error + dev_2 = 0.0 + omega_2 = 0.0 + do alpha = 1,N_monomer + dev_2 = dev_2 + norm_2(dev(itr_prev+1,:,alpha))**2.0 + omega_2 = omega_2 + norm_2(omega(alpha,2:))**2.0 + end do + + error_new_1 = (dev_2/omega_2)**0.5 - itr_stress = itr_stress + 1 - - iterative_loop : do + if(domain)then + dev_2 = dev_2 + norm_2(dev_cp(itr_prev+1,:))**2.0 + omega_2 = omega_2 + norm_2(cell_param(1:N_cell_param))**2.0 + endif + + error_new = (dev_2/omega_2)**0.5 - itr = itr + 1 - if(itr max_itr ) then converge = .false. - write(*,*)'2' exit iterative_loop else if ( error > 100.0*large ) then converge = .false. - write(*,*)'3' exit iterative_loop endif + + if(itr==1)then + do alpha=1,N_monomer + omega(alpha,2:) = omega(alpha,2:) + lam_AM*dev(itr,:,alpha) + end do + if(domain) cell_param(1:N_cell_param) = cell_param(1:N_cell_param) + lam_AM*dev_cp(itr,:) - if(itr==3) write(6,*)'Anderson-Mixing Starts Now ' - - WW = 0.0 - DD = 0.0 - - if(itr 0) then - call mu_phi_solvent(mu_solvent, phi_solvent, q_solvent) - endif - call free_energy(N, rho, omega, phi_chain, mu_chain, & - phi_solvent, mu_solvent, f_Helmholtz, pressure) - !# ifdef DEVEL - call divide_energy(rho,omega,phi_chain,phi_solvent,Q,f_Helmholtz,f_component,overlap) - !# endif - - end do iterative_loop - - - ! AM scheme first resolve omega fields for a fixed unt cell dimenions - ! and then take one step of unit cell dimension. Then again resolves - ! the omega fields for updated unit cell dimensions and conitnue this - ! to iterate fields and unit cell parameterssequentially until both converges. - if(domain) then - - if ( maxval(abs(stress)*stress_rescale) < error_max ) then - converge = .true. - write(*,*)'1' - exit stress_loop - else if ( itr_stress > 20 ) then - converge = .false. - write(*,*)'2' - exit stress_loop - else if ( maxval(abs(stress))*stress_rescale > large ) then - converge = .false. - write(*,*)'3' - exit stress_loop - endif - - write(6,*)'Stress relaxation iteration = ',itr_stress - - call response_dstress_dcell(N,omega,dstress_dcell ) - - if(N_cell_param==1)then - p = -stress(1)/dstress_dcell(1,1) - else - call dgetrf(N_cell_param,N_cell_param,dstress_dcell,N_cell_param,ipvt,info) - if(info/=0) stop "Full dstress_dcell LU factorization failed." - call dgetri(N_cell_param,dstress_dcell,N_cell_param,ipvt,work,lwork,info) - if(info/=0) stop "Full dstress_dcell inversion failed." - - p = - matmul(dstress_dcell,stress) - endif - - vsum=0.0_long - vsum = vsum + dot_product(cell_param(1:N_cell_param),&! - cell_param(1:N_cell_param)) - stpmax = STPMX * max( sqrt(vsum), dble(N_cell_param) ) - - vsum = sqrt(dot_product(p,p)) - if (vsum > stpmax) then - write(6,*) "residual rescaled" - p = p * stpmax / vsum - end if + do alpha=1,N_monomer + elm = 0 + elm = (dev(itr_prev+1,:,alpha) - dev(itr_prev+1-mm,:,alpha))&! + *(dev(itr_prev+1,:,alpha) - dev(itr_prev+1-nn,:,alpha)) + Umn(mm,nn) = Umn(mm,nn) + sum(elm) + end do + if(domain)then + elm_cp = (dev_cp(itr_prev+1,:) - dev_cp(itr_prev+1-mm,:))&! + *(dev_cp(itr_prev+1,:) - dev_cp(itr_prev+1-nn,:)) + Umn(mm,nn) = Umn(mm,nn) + sum(elm_cp) + endif + + Umn(nn,mm) = Umn(mm,nn) + + end do + + do alpha=1,N_monomer + elm = 0 + elm = (dev(itr_prev+1,:,alpha) - dev(itr_prev+1-mm,:,alpha))&! + *dev(itr_prev+1,:,alpha) + Vm(mm,1) = Vm(mm,1) + sum(elm) + end do + if(domain)then + elm_cp = 0 + elm_cp = (dev_cp(itr_prev+1,:) - dev_cp(itr_prev+1-mm,:))&! + *dev_cp(itr_prev+1,:) + Vm(mm,1) = Vm(mm,1) + sum(elm_cp) + endif + + end do + + ! Inverting Umn + call dgetrf(itr_prev,itr_prev,Umn,itr_prev,ipvt,info) + if(info/=0) stop "Full Umn LU factorization failed." + call dgetri(itr_prev,Umn,itr_prev,ipvt,work,lwork,info) + if(info/=0) stop "Full Umn inversion failed." + Cn = matmul(Umn,Vm) + + do alpha = 1,N_monomer + WW = omega(alpha,2:) + DD = dev(itr_prev+1,:,alpha) + do nn=1,itr_prev - do i=1,N_cell_param - cell_param(i) = cell_param(i) + p(i) - end do - - call make_unit_cell - call make_ksq(G_basis) - do i = 1, N_cell_param - call make_dGsq(dGsq(:,i), dGG_basis(:,:,i)) - end do + WW = WW + Cn(nn,1)*(om_hist(itr_prev+1-nn,:,alpha) - om_hist(itr_prev+1,:,alpha)) + DD = DD + Cn(nn,1)*(dev(itr_prev+1-nn,:,alpha) - dev(itr_prev+1,:,alpha)) + + end do + + omega(alpha,2:) = WW + lam_AM*DD + + end do - stress = scf_stress(N, N_cell_param, dGsq) + if(domain)then + WW_cp = cell_param(1:N_cell_param) + DD_cp = dev_cp(itr_prev+1,:) + do nn=1,itr_prev - write(6,*)'Updated Unit Cell dimensions = ',cell_param(1:N_cell_param) - write(6,*)'stress = ',stress - write(6,*)'' - -!! end do stress_relax + WW_cp = WW_cp + Cn(nn,1)*(cp_hist(itr_prev+1-nn,:) - cp_hist(itr_prev+1,:)) + DD_cp = DD_cp + Cn(nn,1)*(dev_cp(itr_prev+1-nn,:) - dev_cp(itr_prev+1,:)) + + end do + + cell_param(1:N_cell_param) = WW_cp + lam_AM*DD_cp - call density(N, omega, rho, q, q_solvent) + endif + + ! Note that Umn, Vm, Cn change sizes after every iteration + ! until N_hist number of iterations are completed, so have + ! to deallocate them everytime and allocate at the start of step + if(itr Date: Sat, 2 Jun 2018 17:34:08 -0500 Subject: [PATCH 8/9] Modified error calculation step in AM to compute error for any number of monomer types (Earlier code was calculating only for systems have 2 and 3 different monomer types --- src/iterate/iterate_mod.fp.f | 28 +++++----------------------- 1 file changed, 5 insertions(+), 23 deletions(-) diff --git a/src/iterate/iterate_mod.fp.f b/src/iterate/iterate_mod.fp.f index de2307a..0fd4d78 100644 --- a/src/iterate/iterate_mod.fp.f +++ b/src/iterate/iterate_mod.fp.f @@ -1108,7 +1108,7 @@ subroutine iterate_AM( & real(long) :: elm(N-1), elm_cp(N_cell_param) real(long) :: dev_2, omega_2 real(long) :: lam_AM, f_old - real(long) :: error1, error2, error_new, error_new_1, error_A, error_B, error_C + real(long) :: error1, error2, error_new, error_new_1 integer :: mm,nn, itr_prev dev=0 @@ -1188,7 +1188,7 @@ subroutine iterate_AM( & if(domain) dev_cp(itr_prev+1,:) = -stress - ! Calculating Error + ! Calculating Error as defined in Matsen AM Papers dev_2 = 0.0 omega_2 = 0.0 do alpha = 1,N_monomer @@ -1205,31 +1205,13 @@ subroutine iterate_AM( & error_new = (dev_2/omega_2)**0.5 - !! Newly added here from git_modified + !! Calculating Error as Max Residuals if(domain)then - if(N_monomer==2)then - error_A = maxval(abs(dev(itr_prev+1,:,1))) - error_B = maxval(abs(dev(itr_prev+1,:,2))) - error1 = max(error_A,error_B) - elseif(N_monomer==3)then - error_A = maxval(abs(dev(itr_prev+1,:,1))) - error_B = maxval(abs(dev(itr_prev+1,:,2))) - error_C = maxval(abs(dev(itr_prev+1,:,3))) - error1 = max(error_A,error_B,error_C) - endif + error1 = maxval(abs(dev(itr_prev+1,:,:))) error2 = maxval(abs(dev_cp(itr_prev+1,:))) error = max(error1, error2*stress_rescale) else - if(N_monomer==2)then - error_A = maxval(abs(dev(itr_prev+1,:,1))) - error_B = maxval(abs(dev(itr_prev+1,:,2))) - error1 = max(error_A,error_B) - elseif(N_monomer==3)then - error_A = maxval(abs(dev(itr_prev+1,:,1))) - error_B = maxval(abs(dev(itr_prev+1,:,2))) - error_C = maxval(abs(dev(itr_prev+1,:,3))) - error1 = max(error_A,error_B,error_C) - endif + error1 = maxval(abs(dev(itr_prev+1,:,:))) error = error1 endif From 7805065796f7351cafdfed2c8cbcc26cd104183b Mon Sep 17 00:00:00 2001 From: David Morse Date: Sun, 3 Jun 2018 12:23:41 -0500 Subject: [PATCH 9/9] Editing user guide. --- doc/user-man/compile-make.rst | 30 +++++++-------- doc/user-man/param.rst | 69 ++++++++++++++++++++--------------- 2 files changed, 54 insertions(+), 45 deletions(-) diff --git a/doc/user-man/compile-make.rst b/doc/user-man/compile-make.rst index c6b4b5a..c92b063 100644 --- a/doc/user-man/compile-make.rst +++ b/doc/user-man/compile-make.rst @@ -5,16 +5,16 @@ Compiling from source, using make ================================= It is also possible to compile using the unix make utility alone, using -a simple Makefile that is provided in the make/ directory of the git +a simple Makefile that is provided in the make/ subdirectory of the git repository. The instructions for using make to compile from source are the same on any unix-like operating system, including Max OS X. The main -difference among different unix environments is the locations of the -required libraries. +differences among different unix environments are the locations of the +required libraries. To compile the code in this way, proceed as follows: * Follow the directions given in the discussion of - :ref:`install-compile-cmake-dependencies-sub` on the previoius + :ref:`install-compile-cmake-dependencies-sub` on the previous page. You will need to install all dependencies listed there except cmake. @@ -22,15 +22,15 @@ To compile the code in this way, proceed as follows: :ref:`install-compile-cmake-getsource-sub` on the previous page to create an appropriate directory structure and obtain the source code. After this step, you should have a directory named - pscf/ with a directory named git/ that contains the contents of - the git repository. You do not need to create a subdirectory - of pscf/ named cmake/ if you are not using cmake. + pscf/ with a subdirectory named git/ that contains the contents + of the git repository. If you are not using cmake, then you do + not need to create a subdirectory of pscf/ named cmake/. - * Change the working directory (cd) to the directory pscf/git/make . + * Change working directory (cd) to the directory pscf/git/make . Note that this is an existing subdirectory of the pscf/git directory, and is different from the initially empty directory - pscf/cmake from which we recommended that invoke cmake when using - cmake to compile. + pscf/cmake from which we recommended cmake be invoked when + using cmake to compile. * The pscf/git/make directory will contain files named config.mk_r and Makefile. Make a copy of the file config.mk_r, by entering:: @@ -49,8 +49,8 @@ To compile the code in this way, proceed as follows: > make -j4 all from within pscf/git/make. The "-j4" option is not necessary, and - simply installs make to try to use up to 4 CPU cores during - compilation, if available. + simply instructs the make utility to try to use up to 4 CPU cores + during compilation, if multiple cores are available. * To install in the directory specified by the $(INSTALL) makefile variable (as defined in config.mk), enter:: @@ -66,9 +66,9 @@ Several of these steps are discussed in more detail below. **Editing the config.mk configfuration file** In the config.mk file in the src/make directory (which you should have -created by copying config.mk_r), you will need to set values for several -macro variables to values appropriate to your system. Makefile variables -you may need to reset are: +created by making a copy of config.mk_r), you will need to set values for +several macro variables to values appropriate to your system. Makefile +variables you may need to reset are: ========= ======================================================== F90 path to executable for Fortran 90 compiler diff --git a/doc/user-man/param.rst b/doc/user-man/param.rst index 21385d7..ac067ec 100644 --- a/doc/user-man/param.rst +++ b/doc/user-man/param.rst @@ -5,16 +5,17 @@ Parameter File ************** -The main program reads an parameter file containing the parameters and -instructions for a calculation. This file is divided into sections, -each of which contains a different type of information. +The main program reads a parameter file containing the parameters and +instructions for a calculation. This file is divided into sections, each +of which contains a different type of information. Each section is preceded by a blank line and starts with a line containing a section title in all capital letters (i.e., 'CHEMISTRY', 'UNIT_CELL', -etc.) Each section may contain values of a sequence of variables. The -name of each variable appears on a line by itself, followed by the value -or (for arrays) values on one more subsequent lines. The variables within -each section must appear in a predetermined (i.e., hard-coded) order. +etc.) Each section may contain values of a sequence of variables. The +variables within each section must appear in a predetermined (i.e., hard-coded) +order. The name of each variable appears on a line by itself, followed by +either the variable value or, for array-valued variables, by a sequence of +values of the elements of the array on one more subsequent lines. The program stops when it encounters the section title 'FINISH'. @@ -156,22 +157,28 @@ parameter file section is given below in a discussion of :ref:`param-sections-se :ref:`param-finish-sub` Stop program =============================== ==================================================== +**Common Workflows** + Several standard types of computation are possible using the blocks listed above: - - Iterate: To solve solve SCF equations for a single state point, include - all of the listed below sections except the SWEEP and RESPONSE sections. + - Iterate: To solve solve the SCF equations for a single state point, include + all of the sections above, in the order listed, except SWEEP and RESPONSE + sections, which should not appear. Also exclude the SOLVENTS section if the + system of interest is a polymer melt or a mixture of polymers with no small + molecule solvent component. - Sweep: To compute a sequence of different states along a line in parameter - space, include both an ITERATE and SWEEP function, but not a RESPONSE - section. The ITERATE section must precede the SWEEP section, and is used - to obtain a solution for the initial choice of parameters. + space, include all of the sections listed above, in the order listed, except + the RESPONSE section. The ITERATE section must precede the SWEEP section, + and is used to obtain a solution for the initial choice of parameters, which + is then used as a starting solution for the rest of the sweep. - - Response: To compute the self-consistent-field or RPA linear susceptibility of a - periodic microstructure, include ITERATE and RESPONSE sections, but do not include - a SWEEP section. + - Response: To compute the self-consistent-field or RPA linear susceptibility + of a periodic microstructure, include ITERATE and RESPONSE sections, but do + not include a SWEEP section. -The SOLVENTS section may be omitted for calculations on polymer melts, with no small -molecule solvent. +The SOLVENTS section may always be omitted for calculations on systems that do +not contain any small molecule solvent. **Miscellaneous Utilities** @@ -182,13 +189,14 @@ transformations on fields or parameters, or to output additional information. Section Description ============================== =============================================== :ref:`param-fieldtorgrid-sub` Read field file in symmetry-adapated format - and output file in coordinate grid format + and output in coordinate grid format :ref:`param-rgridtofield-sub` Read field in coordinate grid file format and output in symmetry-adapated format - :ref:`param-kgridtorgrid-sub` Read field in k-space and output in r-space + :ref:`param-kgridtorgrid-sub` Read field in k-space grid format and output + in r-space coordinate grid format :ref:`param-rhotoomega-sub` Read rho field, compute and output omega field :ref:`param-rescale-sub` Redefine monomer reference volume - :ref:`param-waves-sub` Output map of waves to basis functions + :ref:`param-waves-sub` Output relationship of waves to basis functions :ref:`param-group-sub` Output all elements of space group ============================== =============================================== @@ -226,7 +234,7 @@ of the solvent volume to the reference volume. Note that PSCF does not require the user to input a value for the monomer reference volume - the choice of reference volume is implicit in the values -given for other quantities. Changes in ones choice of reference volume lead +given for other quantities. Changes in one's choice of reference volume lead to corresponding changes in the values for the chi parameters, which are proportional to the reference volume, and in the kuhn lengths, which are proportional to the square root of the reference volume. @@ -588,7 +596,7 @@ If "itr_algo" is "NR", PSCF uses a quasi-Newton-Raphson iteration algorithm that is unique to this program. This algorithm constructs a physically motivated initial approximation for the Jacobian matrix in which elements associated with long wavelength components of the -:\math:`\omega` field are computed numerically and shorter wavelength +:math:`\omega` field are computed numerically and shorter wavelength components are estimated. After construction and inversion of this initial estimate, Broyden updates of the inverse Jacobian are used to refine the estimate of the inverse Jacobian. This method requires a @@ -773,7 +781,7 @@ the input variable "vref_scale". Variable Type Description ================ ============= ============================ input_filename character(60) input omega file name - vref_scale real scale factor + vref_scale real scale factor :math:`\lambda' ================ ============= ============================ This command applies the following set of transformations to each @@ -808,19 +816,20 @@ rescaled choice of parameter values. Using the RESCALE command to read in a file containing such a converged solution should thus cause the subsequent ITERATE command to terminate immediately, since the error should be less than the numerical threshhold on -and output the new parameters to an output summary file and the -rescaled omega field to an output omega file. +input. This would cause the program to immediately output the +rescaled parameters to an output summary file and the rescaled +omega field to an output omega file. .. _param-waves-sub: OUTPUT_WAVES ------------ -Output the relationship between plane waves and symmetry-adapted -basis functions, by outputting a file containing showing which -star each wavevector belongs to and the coefficients of the -plane-wave within a symmetry adapted basis function assocated -with that star. +This command outputs a file that describes the relationship between +complex exponential plane wave basis functions and symmetry-adapted +basis functions. The resulting file lists which star each wavevector +belongs to and the coefficient of the plane-wave within a symmetry +adapted basis function assocated with that star. ================ ============= ============================ Variable Type Description