diff --git a/README.md b/README.md index 05ba3c9ec..f89a39a6d 100644 --- a/README.md +++ b/README.md @@ -29,7 +29,7 @@ functions or small PAO basis sets are also possible. CONQUEST reads pseudopotentials produced by Don Hamann's [ONCVPSP](http://www.mat-simresearch.com) code, and is fully -compatible with the the [PseudoDojo](http://www.pseudo-dojo.org) +compatible with the [PseudoDojo](http://www.pseudo-dojo.org) and [SG15](http://www.quantum-simulation.org/potentials/sg15_oncv/) databases. It can also read the .ion files produced by Siesta, using both the pseudopotentials and @@ -38,7 +38,7 @@ pseudo-atomic orbitals (PAOs) in those files. CONQUEST integrates with The manual is provided with the release (see the ``docs`` directory) but is also available on -[ReadTheDocs](https://conquest.readthedocs.io/en/latest/). +[ReadTheDocs](https://conquest.readthedocs.io/en/latest/). ## DOI @@ -50,4 +50,4 @@ DOI for the specific version you have used in any given study. ## Version number -CONQUEST is now at version 1.5 (tag:v1.5) +CONQUEST is now at version 1.6 (tag:v1.6) diff --git a/docs/ase-conquest.rst b/docs/ase-conquest.rst index 79f5fb51f..bda951d3d 100644 --- a/docs/ase-conquest.rst +++ b/docs/ase-conquest.rst @@ -24,7 +24,7 @@ set generation executable ``MakeIonFiles`` and pseudopotential database. These variables are: * ``ASE_CONQUEST_COMMAND``: the Conquest executable command including MPI/openMPI prefix. -* ``CQ_PP_PATH``: the PAO path directory to where are located the the ``.ion`` files. +* ``CQ_PP_PATH``: the PAO path directory containing the ``.ion`` files. * (optional) ``CQ_GEN_BASIS_CMD`` : the PAO generation executable ``MakeIonFiles``. .. @@ -73,7 +73,7 @@ for example: 'C' : {'file': 'C_SZP.ion'}} In this case they are all assumed to be obtained from Hamann pseudopotentials, -which are the default. Knowing the the exchange and correlation functional +which are the default. Knowing the exchange and correlation functional from the Conquest input (*vide infra*) and the chemical symbol , the *Calculator* will search the ``.ion`` file in different places: :: @@ -167,17 +167,24 @@ or, equivalently, In basic calculate mode (compute ``energy``), the *Calculator* comes with 3 *methods*: -* ``write_input()``: - this function will setup the :ref:`input files `. For CONQUEST, the PAO basis will be generated/copied with respect to the dictionary key/value pairs, and ``Conquest_input`` file including the calculation parameters will be written, a long with the ``coordinate`` file, containing the lattice vectors (in Bohr Unit) and atomic positions (in fractional coordinates). +* ``write_input()``: + This function sets up the :ref:`input files `. For CONQUEST, the + PAO basis is generated or copied according to the dictionary key/value pairs. + The ``Conquest_input`` file containing the calculation parameters is written + along with the ``coordinate`` file, containing the lattice vectors (in bohr) and + atomic positions (in fractional coordinates). -* ``execute()``: - this function execute the calculation. For CONQUEST, it will launch the ``ASE_CONQUEST_COMMAND`` setup in :ref:`the environment varaibles `. +* ``execute()``: + This function executes the calculation. For CONQUEST, it launches the + ``ASE_CONQUEST_COMMAND`` set in :ref:`the environment variables `. -* ``read_results()``: - this function post-process the the output file. For CONQUEST, the ``energy``, ``forces``, ``stress`` and ``eigenvalues`` will be extracted from the ``Conquest_out_ase`` output file. +* ``read_results()``: + This function post-processes the output file. For CONQUEST, the ``energy``, + ``forces``, ``stress`` and ``eigenvalues`` are extracted from the + ``Conquest_out_ase`` output file. .. note:: - The funtion ``read_results()`` operate on the ``Conquest_out_ase`` file. + The function ``read_results()`` operates on the ``Conquest_out_ase`` file. This output file is not created by default by CONQUEST. If you want to post-process a calculation with an input generated by hand you must add ``IO.WriteOutToASEFile True`` in ``conquest_input``. @@ -266,7 +273,7 @@ keyword type default value description ``IO.WriteOutToASEFile`` ``bool`` True write ASE output file ; **must always be True when using ASE for post-processing** ``IO.Iprint`` ``int`` 1 verbose for the output ; **must always be 1 when using ASE for post-processing** ``DM.SolutionMethod`` ``str`` 'diagon' 'diagon' stands for diagonalisation other is 'ordern' (base on density matrix) -``General.PseudopotentialType`` ``str`` 'Hamann' kind of pseudopotential other type are 'siesta' and 'abinit' +``General.PseudopotentialType`` ``str`` 'Hamann' pseudopotential type; accepted values are 'Hamann', 'plato' and 'siesta' ``SC.MaxIters`` ``int`` 50 maximum number SCF cycles ``AtomMove.TypeOfRun`` ``str`` 'static' 'static' stands for single (non)SCF other are 'md' or optimisation algorithms. ``Diag.SmearingType`` ``int`` 1 1 for Methfessel-Paxton ; 0 for Fermi-Dirac @@ -452,14 +459,15 @@ the first K/L matrix, then adding the following keywords to the calculator: :: - other_keywords = {"General.LoadL": True, + other_keywords = {"General.LoadDM": True, "SC.MakeInitialChargeFromK": True} -These keywords respectively cause the K or L matrix to be loaded from file(s) -``Kmatrix.i**.p*****``, and the initial charge density to be constructed from -this matrix. In all subsequent calculations, the K or L matrix will be written -at the end of the calculation and used as the initial guess for the subsequent -ionic step. +These keywords respectively cause the K or L matrix to be loaded from files +``Kmatrix2.i**.p******`` or ``Lmatrix2.i**.p******``, together with the +corresponding ``InfoGlobal.i**`` metadata file, and the initial charge density +to be constructed from this matrix. In all subsequent calculations, the K or L +matrix will be written at the end of the calculation and used as the initial +guess for the subsequent ionic step. Go to :ref:`top `. @@ -514,5 +522,3 @@ interpolates the equation of state and uses ``matplotlib`` to generate a plot. eos.plot('diamond-eos.pdf') # Plot the equation of state Go to :ref:`top `. - - diff --git a/docs/conf.py b/docs/conf.py index 55aaf60f9..2f99118c1 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -63,9 +63,9 @@ # built documents. # # The short X.Y version. -version = u'1.5' +version = u'1.6' # The full version, including alpha/beta/rc tags. -release = u'1.5' +release = u'1.6' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/docs/errors.rst b/docs/errors.rst index 488dc94fa..84b818daa 100644 --- a/docs/errors.rst +++ b/docs/errors.rst @@ -1,10 +1,258 @@ .. _errors: -=========== -Error codes -=========== - -Error codes generated by different CONQUEST routines will be collated -here, along with explanations of the root cause and suggested fixes. -This will form part of the full release, but is not implemented in the -pre-release. +========================== +Errors and troubleshooting +========================== + +CONQUEST reports most fatal errors as text rather than as numbered error +codes. This page lists common messages for which the immediate cause and a +safe first check can be identified. It is not a complete catalogue of every +internal consistency check in the code. + +Where a message contains ````, ```` or ````, the angle- +bracketed text stands for a value supplied at runtime; the surrounding text +matches the source diagnostic. + +How CONQUEST reports problems +----------------------------- + +Fatal errors are written to ``stdout`` and to the configured main output unit +before CONQUEST stops the MPI calculation. The main output is normally +``Conquest_out``, but its destination can be changed; see :ref:`io_output_main`. +For an early startup failure, the output file might not yet have been opened. +Do not keep only the final MPI abort line: the CONQUEST message immediately +before it is normally more useful. + +Warnings issued through CONQUEST's warning routine are written to +``Conquest_warnings`` and, depending on the ``IO.Iprint`` level, can also +appear in the main output. Check ``Conquest_warnings`` after every run, +including runs that appear to finish normally. + +A ``CQ.stop`` file requests that a running calculation stop when CONQUEST +next checks for it. A stop caused by that file is deliberate and should not +be diagnosed as an error. + +First checks after a failed run +------------------------------- + +#. Preserve the complete error text, ``Conquest_out``, + ``Conquest_warnings`` and any scheduler output before restarting. +#. Find the first CONQUEST error or warning, rather than relying only on the + last message printed by MPI or the scheduler. +#. Check that the run directory contains the intended input, coordinate, + ion/pseudopotential and restart files. +#. Check coordinate conventions, distance units and species numbers against + :ref:`io_coords` and :ref:`input_atomic_spec`. +#. Identify whether the problem is an input error, a numerical convergence + failure, or an external runtime failure before changing tolerances, + iteration limits or the parallel layout. + +Startup and input files +----------------------- + +``We need Conquest_input to run !`` + CONQUEST could not open ``Conquest_input`` in the run directory. Check the + filename, working directory and read permissions. The required input files + and the main input structure are described in :ref:`io_files`. + +``No coordinate file specified: please set with IO.Coordinates`` + No coordinate filename was supplied. Set ``IO.Coordinates`` to the + intended file. A PDB input that cannot be opened instead gives + ``Reading pdb file: file error``. Check the coordinate filename and the + selected format against :ref:`io_coords` and :ref:`input_coords`. + +``Too few species in ChemicalSpeciesLabel:`` + The ``ChemicalSpeciesLabel`` block contains fewer entries than + ``General.NumberOfSpecies``. Make the count and block agree and ensure that + the block is terminated correctly; see :ref:`input_atomic_spec`. + +``Species specified greater than number in input file:`` + A species index in the coordinate file is greater than the number of + species declared in the input. Check the coordinate indices and the + numbered entries in ``ChemicalSpeciesLabel``. See + :ref:`input_atomic_spec` for the block format. + +``read_pseudopotential: file error`` + A file used by the legacy pseudopotential reader could not be opened. + Check the filename printed immediately before this error, its location and + its permissions. Do not substitute an unrelated pseudopotential merely to + pass this check; see :ref:`io_ion`, :ref:`input_general` and :ref:`basissets`. + +``Functionals differ between pseudopotential files:`` + When ``General.FunctionalType`` is left at its default of zero, CONQUEST + derives the functional from the first pseudopotential and finds a different + one in another file. Use a mutually consistent set of files; see + :ref:`input_general` and :ref:`io_ion`. + +``Functional in input file differs to pseudopotential:`` + ``General.FunctionalType`` differs from the functional stored in a + pseudopotential. The calculation stops unless + ``General.DifferentFunctional`` explicitly permits the mismatch. Prefer + matching the input and pseudopotentials; an override does not make an + incompatible pseudopotential valid. See :ref:`input_general`. + +``Functional in input file differs to pseudopotential but proceeding:`` + The same mismatch was found, but ``General.DifferentFunctional`` permits + the run to continue with a warning. Verify that the override is intended + and that the input and pseudopotentials are physically compatible; see + :ref:`input_general` and :ref:`io_ion`. + +Coordinates and parallel layout +------------------------------- + +``Expected fractional coordinates but many are greater than one:`` + This warning means that several values look Cartesian although + ``IO.FractionalAtomicCoords T`` selects fractional coordinates. + +``Expected Cartesian coordinates but many are less than one:`` + This warning is the converse heuristic: several values look fractional + although ``IO.FractionalAtomicCoords F`` selects Cartesian coordinates. In + either case, verify ``IO.FractionalAtomicCoords``, + ``General.DistanceUnits`` and the coordinate file itself. The warning is a + heuristic, so inspect the data rather than changing the flag automatically; + see :ref:`io_coords` and :ref:`input_coords`. + +``Atoms are too close to each other:`` + At least two non-ghost atoms are at or below the configured minimum initial + separation. + +``Atoms are too far apart! Minimum distance is`` + Even the smallest detected separation is above the configured maximum. + For either distance error, first inspect duplicated atoms, cell dimensions, + coordinate wrapping, ``IO.FractionalAtomicCoords`` and + ``General.DistanceUnits``. Do not disable the check until the geometry and + units have been verified; see :ref:`io_coords` and :ref:`input_coords`. + +``We must have at least one atom per process:`` + The run uses more MPI processes than atoms. Rerun with no more processes + than atoms; see :ref:`io_coords` for the coordinate file's atom count. + +``More processors than partitions!`` + A manual ``General.NPartitionsX/Y/Z`` product is smaller than the MPI + process count. The related messages ``sfc_partitions_module: too few + partitions created for number of processors`` and + ``sfc_partitions_module: too few occupied partitions for number of + processors`` mean that the generated or occupied partition count is also + insufficient. Reduce the process count or review the partition settings + in :ref:`input_general`; do not assume that adding more processes will + accelerate a small or sparse system. + +Diagonalisation and convergence +------------------------------- + +``Diag: proc grid product is zero:`` + A manually specified diagonalisation processor grid has a zero dimension. + Check ``Diag.ProcRows`` and ``Diag.ProcCols`` against :ref:`input_diag`. + +``Diag: proc grid product is too large:`` + ``Diag.ProcRows`` multiplied by ``Diag.ProcCols`` exceeds the processes in + the applicable k-point process group. Use a valid grid or let CONQUEST + choose it. See :ref:`gs_diag_para` and :ref:`input_diag`. + +**Messages:** ``block_size_r not a factor of matrix size !``, +``block_size_c not a factor of matrix size !`` or +``Can't find a good block size: please set manually`` + +These errors concern the ScaLAPACK matrix distribution when matrix padding is +disabled. Check ``Diag.PaddingHmatrix``, ``Diag.BlockSizeR`` and +``Diag.BlockSizeC`` against the restrictions in :ref:`gs_pad` and +:ref:`input_diag`. Change them as a consistent layout, not as independent +tuning values. + +``Code compiled without ELPA! Set Diag.UseELPA F`` + The input requests ELPA but the executable was built without ELPA support. + Use ScaLAPACK by disabling ``Diag.UseELPA``, or use an executable built and + linked with ELPA. See :ref:`gs_diag_elpa` and :ref:`install_compile`. + +``FindEvals: pzhegvx failed for mode with INFO=`` + The selected generalized eigensolver returned a status that CONQUEST treats + as fatal. Despite the ``pzhegvx`` wording, this message can also follow an + ELPA call. Preserve the mode, ``INFO`` value and ``Diag.UseELPA`` setting, + together with the parallel layout and full output. Do not infer a physical + cause from the number alone. Check the diagonalisation setup in + :ref:`gs_diag` before reporting a reproducible failure. + +``ScaLAPACK pzhegvx warning, info=`` + The selected eigensolver returned ``INFO=2`` or ``INFO=4``. CONQUEST + permits these statuses to continue and records a warning. This wording is + also used when ``Diag.UseELPA`` is enabled, so interpret the status using + the selected solver. Check the resulting eigenvalues, occupations and + convergence, and retain the exact status if reporting unexpected results; + see :ref:`gs_diag`. + +``SelfCon/PulayMixSC_spin: too many SCF iterations:`` + The self-consistency cycle reached ``SC.MaxIters`` without satisfying its + convergence test, unless ``SC.ContinueOnSCFail`` was selected. Inspect the + residual history and verify that the basis, grid and physical setup are + sensible before increasing the limit. Mixing and convergence controls are + discussed in :ref:`conv_scf` and :ref:`gs_scf`. + +Atomic movement and restart files +--------------------------------- + +**Messages:** ``Step too small: safemin2 failed!``, +``Step too small: safemin_cell failed!`` or +``Step too small: safemin_full failed!`` + +An ionic, cell or combined line search reduced its trial step below the usable +threshold. Check the starting geometry, forces, preceding energy changes and +line-search messages before changing optimisation controls. See +:ref:`strucrelax` for the supported relaxation methods and their inputs. + +**Messages:** ``Error: Too many SHAKE iterations !`` or +``Error: Too many RATTLE iterations !`` + +A rigid-bond position or velocity correction did not converge within its +iteration limit. Check the constraint definitions, the initial constrained +bond lengths and the preceding dynamics before changing tolerances or limits. +See :ref:`moldyn` for molecular-dynamics setup. + +``Fail in opening InfoGlobal.dat`` + The restart reader could not open the indexed ``InfoGlobal`` metadata file. + Despite the diagnostic's ``.dat`` wording, check the actual filename printed + immediately before it (such as ``InfoGlobal.i00``), the restart settings and + read permissions. See :ref:`md_restart`. + +**Messages:** ``Fail in opening Lmatrix file in Binary.`` or +``Fail in opening ASCII Lmatrix file.`` + +The matrix restart reader could not open a file in the selected binary or text +format. The ``Lmatrix`` wording is also used when reading other matrices, +including ``Kmatrix2``. Check that the complete indexed matrix-file set is +present and that ``IO.MatrixFile.BinaryFormat.Grab`` matches its format. Do +not assume the MPI process count must match the previous run: this reader uses +the process count stored in ``InfoGlobal``. See :ref:`gs_scf_restart` and +:ref:`md_restart`. + +``grab_blip_coeffs: failed to open input file `` + With a blip basis and ``Basis.LoadCoeffs T``, a per-process + ``blip_coeffs`` file could not be opened. Check the filename printed by + CONQUEST, the coefficient files from the same calculation and their read + permissions; see :ref:`basis_readcoeffs` and :ref:`md_restart`. + +External runtime failures +------------------------- + +MPI launch failures, scheduler cancellation, out-of-memory termination, +filesystem errors and Unix signals are not CONQUEST error codes, and their +meaning depends on the platform. Check the scheduler and system logs as well +as the CONQUEST output. For a segmentation fault seen only with multiple +OpenMP threads, also check the thread-stack guidance in :ref:`install_compile`. + +Reporting a reproducible problem +-------------------------------- + +If the checks above do not identify the problem, follow :ref:`faq_bug` and +include the complete diagnostic and enough material to reproduce the run: + +* ``Conquest_input``, the coordinate file and all applicable ion or + pseudopotential files; +* any restart files needed to reproduce a restart failure; +* ``Conquest_out``, ``Conquest_warnings`` and scheduler output; +* the CONQUEST release, commit or other precise version identifier; +* compiler, numerical-library and MPI implementation versions; and +* MPI process and OpenMP thread counts. + +Remove confidential data if necessary, but do not omit the lines immediately +before the failure or silently replace the failing inputs with different +ones. diff --git a/docs/examples.rst b/docs/examples.rst index 0643caf40..00f4857ba 100644 --- a/docs/examples.rst +++ b/docs/examples.rst @@ -277,9 +277,6 @@ We recommend that you work through, in order, the tutorials included in the distribution in the ``tutorials/`` directory to become familiar with the modes of operation of the code. -**NOTE** In the initial pre-release of CONQUEST (January 2020) we have -not included the tutorials; they will be added over the coming months. - Go to :ref:`top `. .. _ex_next: diff --git a/docs/ext-tools.rst b/docs/ext-tools.rst index 887352cbc..ea25a39b9 100644 --- a/docs/ext-tools.rst +++ b/docs/ext-tools.rst @@ -15,176 +15,91 @@ densities, DOS and STM images in useful forms. It is described fully :ref:`here .. _et_md_scripts: -Molecular dynamics analysis +Molecular-dynamics analysis --------------------------- -Several scripts that may be helpful with postprocessing molecular dynamics are -included with CONQUEST. The can be found in the ``tools`` directory, and the -executables are ``plot_stats.py``, ``md_analysis.py`` and ``heat_flux.py``. They -have the following dependencies: +The consolidated ``md_analysis.py`` utility is located in +``src/utilities``. It can plot the thermodynamic statistics written by +CONQUEST and analyse trajectories for radial distribution functions (RDFs), +velocity autocorrelation functions (VACFs), mean-squared displacements (MSDs) +and stress. It can also compare the statistics from several calculations and +calculate heat-flux autocorrelation functions (HFACFs). - * Python 3 - * Scipy/Numpy - * Matplotlib +The utility requires Python 3, NumPy, SciPy and Matplotlib. Run it from the +calculation directory, giving the path to the copy in the CONQUEST source tree, +for example: -If Python 3 is installed the modules can be added easily using ``pip3 install -scipy`` etc. - -These scripts should be run in the calculation directory, and will automatically -parse the necessary files, namely ``Conquest_input``, ``input.log``, -``md.stats`` and ``md.frames`` assuming they have the default names. They will -also read the CONQUEST input flags to determine, for example, what ensemble is -used, and process the results accordingly. +:: -Go to :ref:`top `. + python3 /path/to/CONQUEST/src/utilities/md_analysis.py [options] -.. _et_plot_stat: +The utility reads ``Conquest_input`` and the coordinate file named there. By +default, the MD data are read from ``md.stats``, ``md.frames`` and, for an +HFACF, ``md.heatflux``. A different statistics or frames file can be selected +with ``--stats FILE`` or ``--frames FILE``. CONQUEST controls how often +frames are written to ``md.frames`` with ``AtomMove.OutputFreq``; statistics +are written to ``md.stats`` at every MD step. Plotting statistics +++++++++++++++++++ -:: - - usage: plot_stats.py [-h] [-c] [-d DIRS [DIRS ...]] - [--description DESC [DESC ...]] [--skip NSKIP] - [--stop NSTOP] [--equil NEQUIL] [--landscape] - [--mser MSER_VAR] - - Plot statistics for a CONQUEST MD trajectory - - optional arguments: - -h, --help show this help message and exit - -c, --compare Compare statistics of trajectories in directories - specified by -d (default: False) - -d DIRS [DIRS ...], --dirs DIRS [DIRS ...] - Directories to compare (default: .) - --description DESC [DESC ...] - Description of graph for legend (only if using - --compare) (default: ) - --skip NSKIP Number of equilibration steps to skip (default: 0) - --stop NSTOP Number of last frame in analysis (default: -1) - --equil NEQUIL Number of equilibration steps (default: 0) - --landscape Generate plot with landscape orientation (default: - False) - --mser MSER_VAR Compute MSER for the given property (default: None) - -Running ``plot_stats.py --skip 200`` in your calculation will generate a plot -which should resemble the example below, skipping the first 200 steps. This -example is a molecular dynamics simulation of 1000 atoms of bulk silicon in the -NPT ensemble, at 300 K and 0.1 GPa. +With no analysis option, the utility plots the energy contributions, conserved +quantity, temperature and pressure from ``md.stats`` and writes ``stats.pdf``. +For an NPT calculation it also plots the volume. For example, -.. image:: stats.jpg +:: -The four plots are respectively the breakdown of energy contributions, the -conserved quantity, the temperature and the pressure, the last of which is only -included for NPT molecular dynamics. Several calculations in different -directories can be compared using ``plot_stats.py --compare -d dir1 -dir2 --description "dir1 description" "dir2 description"``. The following -example compares the effect of changing the L tolerance in the above simulation. -Note that the contents of the description field will be in the legend of the -plot. + python3 /path/to/CONQUEST/src/utilities/md_analysis.py --skip 200 -.. image:: compare.jpg - -Go to :ref:`top `. +omits the first 200 MD steps from the plots. ``--equil N`` independently sets +the initial part omitted when calculating averages, while ``--landscape`` uses +a two-by-two plot layout. -.. _et_md_ana: +.. image:: stats.jpg -MD analysis -+++++++++++ +Statistics from calculations in several directories can be compared with +``--compare``. Supply the directories with ``--dirs`` and, optionally, one +legend label per directory with ``--description``: :: - usage: md_analysis.py [-h] [-d DIRS [DIRS ...]] [--skip NSKIP] - [--stride STRIDE] [--snap SNAP] [--stop NSTOP] - [--equil NEQUIL] [--vacf] [--msd] [--rdf] [--stress] - [--nbins NBINS] [--rdfwidth RDFWIDTH] [--rdfcut RDFCUT] - [--window WINDOW] [--fitstart FITSTART] [--dump] - - Analyse a CONQUEST MD trajectory - - optional arguments: - -h, --help show this help message and exit - -d DIRS [DIRS ...], --dirs DIRS [DIRS ...] - Directories to compare (default: .) - --skip NSKIP Number of equilibration steps to skip (default: 0) - --stride STRIDE Only analyse every nth step of frames file (default: - 1) - --snap SNAP Analyse Frame of a single snapshot (default: -1) - --stop NSTOP Number of last frame in analysis (default: -1) - --equil NEQUIL Number of equilibration steps (default: 0) - --vacf Plot velocity autocorrelation function (default: - False) - --msd Plot mean squared deviation (default: False) - --rdf Plot radial distribution function (default: False) - --stress Plot stress (default: False) - --nbins NBINS Number of histogram bins (default: 100) - --rdfwidth RDFWIDTH RDF histogram bin width (A) (default: 0.05) - --rdfcut RDFCUT Distance cutoff for RDF in Angstrom (default: 8.0) - --window WINDOW Window for autocorrelation functions in fs (default: - 1000.0) - --fitstart FITSTART Start time for curve fit (default: -1.0) - --dump Dump secondary data used to generate plots (default: - False) - -The script ``md_analysis.py`` script performs various analyses of the trajectory -by parsing the `md.frames`` file. So far, these include the radial distribution -function, the velocity autocorrelation function, the mean squared deviation, and -plotting the stress. For example, the command, - -``md_analysis.py --rdf --stride 20 --rdfcut 8.0 --nbins 100 --dump --skip 200 --stop 400`` - -computes the radial distribution function of the simulation in the first example -from every 20th time step (every 10 fs in this case), stopping after 400 steps, -with a cutoff of 8.0 A, and the histogram is divided into 100 bins. + python3 /path/to/CONQUEST/src/utilities/md_analysis.py \ + --compare --dirs dir1 dir2 \ + --description "first calculation" "second calculation" -.. image:: rdf.jpg +.. image:: compare.jpg -Go to :ref:`top `. +Trajectory analysis ++++++++++++++++++++ -.. _et_cq_struc: +The options ``--rdf``, ``--vacf``, ``--msd`` and ``--stress`` select analyses +of ``md.frames``. ``--stride N`` uses every Nth frame, ``--skip N`` ignores +the initial part of the trajectory and ``--stop N`` sets the final MD step to +analyse. ``--dump`` writes the numerical data used for RDF, VACF and MSD plots +in addition to the PDF output. -CONQUEST structure file analysis -++++++++++++++++++++++++++++++++ +For example, :: - usage: structure.py [-h] [-i INFILE] [--bonds] [--density] [--nbins NBINS] - [-c CUTOFF [CUTOFF ...]] [--printall] + python3 /path/to/CONQUEST/src/utilities/md_analysis.py \ + --rdf --stride 20 --rdfcut 8.0 --rdfwidth 0.08 \ + --dump --skip 200 --stop 400 - Analyse a CONQUEST-formatted structure +calculates an RDF with an 8 Angstrom cutoff and 0.08 Angstrom bin width, using +every twentieth selected frame. The results are written to ``rdf.pdf`` and, +because ``--dump`` is present, ``rdf.dat``. - optional arguments: - -h, --help show this help message and exit - -i INFILE, --infile INFILE - CONQUEST format structure file (default: - coord_next.dat) - --bonds Compute average and minimum bond lengths (default: - False) - --density Compute density (default: False) - --nbins NBINS Number of histogram bins (default: 100) - -c CUTOFF [CUTOFF ...], --cutoff CUTOFF [CUTOFF ...] - Bond length cutoff matrix (upper triangular part, in - rows (default: None) - --printall Print all bond lengths (default: False) - -The script ``structure.py`` can be used to analyse a CONQUEST-formatted -structure file. This is useful to sanity-check the bond lengths or density, -since an unphysical structure is so often the cause of a crash. For example, the -bond lengths can be computed with - -``structure.py --bonds -c 2.0 3.0 3.0`` - -where the ``-c`` flag specifies the bond cutoffs for the bonds 1-1, 1-2 and 2-2, -where 1 is species 1 as specified in ``Conquest_input`` and 2 is species 2. The -output will look something like this: +.. image:: rdf.jpg -:: +The VACF and MSD analyses similarly write ``vacf.pdf`` and ``msd.pdf``, with +``vacf.dat`` and ``msd.dat`` requested by ``--dump``. ``--stress`` writes +``stress.pdf``; for an NPT trajectory the cell dimensions are plotted with the +stress. - Mean bond lengths: - O-Si: 1.6535 +/- 0.0041 (24 bonds) - Minimum bond lengths: - O-Si: 1.6493 +For heat-flux analysis, enable ``MD.HeatFlux`` in the CONQUEST calculation and +use ``--hfacf --acfwindow T``, where ``T`` is the autocorrelation-window length +in femtoseconds. The utility reads ``md.heatflux`` and writes ``hfacf.pdf``. Go to :ref:`top `. @@ -199,7 +114,7 @@ atomistic simulations. ASE contains a CONQUEST interface, also called *Calculator* so that it can be used to calculate ``energies``, ``forces`` and ``stresses`` as inputs to other calculations such as `Phonon `_ or `NEB `_ that -are not implemented in CONQUEST. ASE is a versatil tool to manage CONQUEST +are not implemented in CONQUEST. ASE is a versatile tool to manage CONQUEST calculations without pain either: * in a **direct** way where pre-processing, calculation and post-processing are managed on-the-fly by ASE, diff --git a/docs/groundstate.rst b/docs/groundstate.rst index 360739535..77dc8ad88 100644 --- a/docs/groundstate.rst +++ b/docs/groundstate.rst @@ -156,10 +156,13 @@ finite gap, you can force integer occupancies with the flag: This is false by default, and will lead to significant problems if the system is metallic at any point during the self-consistent ground state search, so should be considered carefully. -In both of these cases, if there is gap in the band structure, it is written out along with -the location in the Brillouin zone where it is found (though this result will depend very -sensitively on the k-point sampling used). The Fermi level will be placed at the mid-point -between the valence band maximum and the conduction band minimum. +With Fermi-Dirac occupations, if a gap is detected it is written out along with +the location in the Brillouin zone where it is found (though this result will +depend very sensitively on the k-point sampling used). By default, the Fermi +level is then moved to the midpoint between the valence-band maximum and the +conduction-band minimum. This adjustment can be disabled with +``Diag.AdjustEf F``. Integer occupations always use the midgap Fermi level, +regardless of ``Diag.AdjustEf``. The Methfessel-Paxton approach :cite:`g-Methfessel:1989ny` to occupations allows much higher smearing temperatures with minimal effect on the free energy (and @@ -423,13 +426,16 @@ The relevant input tags are: 1 3 2 1 0.0772 %endblock -The ``DFTplusU`` block specifies the projectors onto the correlated subspace (at present -these must be PAOs) and requires one line per species. The block specifies n, l, zeta -and the value of U (in Hartrees). We strongly recommend that PAOs created with -equal energies are used (in the PAO generation, add ``Atom.Cutoffs energy`` to the -atom block), and that the largest radius zeta function (i.e. zeta=1) is used. Careful -testing is needed for the projectors and an appropriate value of U (we note that when using -the same pseudopotentials, CONQUEST and QuantumEspresso calculations agree very well). +The ``DFTplusU`` block specifies the PAO projectors onto the correlated +subspace and requires one line for each species to which a Hubbard correction +is applied. Each line gives the species number, :math:`n`, :math:`l`, zeta +index and the value of U (in Hartree), in that order. DFT+U is not available +with a blip basis. We strongly recommend that PAOs created with equal energies +are used (in the PAO generation, add ``Atom.Cutoffs energy`` to the atom block), +and that the largest-radius zeta function (i.e. zeta=1) is used. Careful +testing is needed for the projectors and an appropriate value of U (we note +that when using the same pseudopotentials, CONQUEST and QuantumEspresso +calculations agree very well). The exchange-correlation (XC) functional is (almost always) set by the ion files used to specify the elements (and generated by the :ref:`PAO generation ` tool). diff --git a/docs/input-output.rst b/docs/input-output.rst index c8ea2427a..6fb149d20 100644 --- a/docs/input-output.rst +++ b/docs/input-output.rst @@ -86,8 +86,7 @@ ion file generation are: ``large``. These are all included in the default input files. -Further fine-grained control can be applied to the basis functions; -this will be documented after the pre-release of CONQUEST. +Further fine-grained control can be applied to the basis functions. .. _ONCVPSP: http://http://www.mat-simresearch.com .. _PseudoDojo: https://www.pseudo-dojo.org/ @@ -167,11 +166,12 @@ the key output flag is given. Further output flags are described in :ref:`input * Charge density * Band-resolved charge density (``IO.outputWF``) - * Density of states (``IO.writeDOS``) + * Density of states, obtained by post-processing ``eigenvalues.dat`` * Atom-projected density of states (``IO.write_proj_DOS``) * Atomic charges, using the Mulliken approach (``IO.AtomChargeOutput``) -The Kohn-Sham eigenvalues are output in the ``eigenvalues.dat`` file. +For diagonalisation calculations, the Kohn-Sham eigenvalues and their +occupancies are output in the ``eigenvalues.dat`` file. The charge densities need post-processing to convert from the standard output format to a file compatible with visualisation (current supported formats include Gaussian CUBE file and OpenDX @@ -214,7 +214,7 @@ files: parameters, atomic positions, velocities, forces, stress). * ``md.checkpoint`` --- data required for MD restart, namely atomic velocities and extended system variables. - * ``md.positions`` --- Atomic coordinates saved at the moment of checkpointing + * ``md.position`` --- atomic coordinates saved at the moment of checkpointing * ``trajectory.xsf`` --- atomic coordinates save in .xsf format, which can be visualised using (for example) VMD, if ``AtomMove.WriteXSF`` is true.. diff --git a/docs/input_tags.rst b/docs/input_tags.rst index 074050995..b46b37bcc 100644 --- a/docs/input_tags.rst +++ b/docs/input_tags.rst @@ -7,7 +7,8 @@ Input tags We have broken down the input tags based on the areas of the code where they apply. For each tag, a default is given. Types of value are specified as: *integer*; -*real*; *boolean*; or *string* (optA/optB are given for string options). +*real*; *boolean*; *string* (optA/optB are given for string options); or +*presence flag* (enabled whenever the tag occurs, regardless of its value). .. contents:: Areas :depth: 1 @@ -21,13 +22,14 @@ General General.NumberOfSpecies (*integer*) Number of species in cell - *default*: none + *default*: 1 -General.PseudopotentialType (*string*) siesta/hamann +General.PseudopotentialType (*string*) Type of pseudopotential (in practice, this defines how the local - part of the pseudopotential is handled) + part of the pseudopotential is handled). The normal values are + ``hamann`` for ONCVPSP-derived data and ``siesta`` for SIESTA ion files. - *default*: hamann (read from ion file) + *default*: hamann General.NeutralAtom (*boolean*) Use neutral atom potential or not (removes need for Ewald sum) @@ -40,7 +42,10 @@ General.FunctionalType (*integer*) parameterisations of the LDA available, as well as three variants of the PBE GGA functional, with numbers given below. - *default*: read from ion file (same as pseudopotentials) + *default*: 0 (take the functional from the ion files) + + If a nonzero value differs from the functional stored in the ion files, + ``General.DifferentFunctional`` must be true or CONQUEST will stop. ========================================= ======= ======================= Functional Keyword Ref @@ -72,7 +77,7 @@ General.EnergyUnits (*string*) Ha/Ry/eV *default*: Ha General.DistanceUnits (*string*) a0/bohr/A - **Output only** Chooses units for distance (Bohr: a0/bohr or Ångströms: A + **Output only** Chooses units for distance (Bohr: a0/bohr or Ångströms: A). *default*: a0 @@ -94,18 +99,13 @@ General.PartitionMethod (*string*) File/Hilbert guaranteed - File — Reads a file (NOT recommended) -General.LoadBalance (*string*) partitions/atoms +General.LoadBalance (*string*) atoms/partitions/supportfunctions Applies to Hilbert above; chooses whether to distribute atoms or partitions - evenly between processors (you are *strongly* recommended to use atoms) + evenly between processors, or to balance the number of support functions + (you are *strongly* recommended to use atoms) *default*: atoms -General.ManyProcessors (*boolean*) - Applies to Hilbert above; chooses method for parallelising Hilbert curve work; - “many” processors here probably means more than two - - *default*: T - General.MaxAtomsPartition (*integer*) Applies to Hilbert above; specifies maximum number of atoms allowed in a partition; triggers extra level of recursion in @@ -128,7 +128,7 @@ General.LoadDM (*boolean*) Specifies whether to load a previous density matrix (K or L depending on whether diagonalisation or linear scaling are selected) from files - *default*: F + *default*: T when ``AtomMove.RestartRun`` is true; F otherwise General.LoadRho (*boolean*) Specifies whether to load a previous charge density from files @@ -215,7 +215,7 @@ Atom.LFDRange (*real*) overlap matrix elements from all atoms within this range will be included in the cluster diagonalisation) - *default*: 0.0 + *default*: value of ``Atom.MultisiteRange`` Go to :ref:`top `. @@ -270,6 +270,61 @@ IO.Dump[Har|XC|PS|ES|Tot]Pot (*boolean*) with XC, PS, ES and Tot) *default*: F + +IO.outputWF (*boolean*) + Writes wavefunction coefficients for post-processing. This output is + available only for static diagonalisation calculations; CONQUEST disables + it with a warning for O(N) or non-static calculations. Bands must be + selected either with ``IO.maxnoWF`` and a ``WaveFunctionsOut`` block or + with the ``IO.min_wf_E`` and ``IO.max_wf_E`` energy limits. + + *default*: F + +IO.maxnoWF (*integer*) + Number of explicitly selected bands whose wavefunction coefficients are + written. A positive value requires a ``WaveFunctionsOut`` block and takes + precedence over energy-range selection. + + *default*: 0 + +WaveFunctionsOut (*block*) + Lists the band numbers to output when ``IO.maxnoWF`` is positive. The + block must contain at least ``IO.maxnoWF`` entries, one band number per + entry. + + *default*: none; required when ``IO.maxnoWF`` is positive + +IO.min_wf_E (*real*) + Lower energy limit, in Hartree, for wavefunction output when no explicit + band list is active. For ordinary ``IO.outputWF`` output, the default is + zero. When ``IO.write_proj_DOS`` is true, an omitted limit instead extends + to the bottom of the eigenspectrum. + + *default*: 0, or the bottom of the eigenspectrum for projected-DOS output + +IO.max_wf_E (*real*) + Upper energy limit, in Hartree, for wavefunction output when no explicit + band list is active. For ordinary ``IO.outputWF`` output, the default is + zero. When ``IO.write_proj_DOS`` is true, an omitted limit instead extends + to the top of the eigenspectrum. + + *default*: 0, or the top of the eigenspectrum for projected-DOS output + +IO.WFRangeRelative (*boolean*) + Interprets ``IO.min_wf_E`` and ``IO.max_wf_E`` relative to the Fermi level + when true, and as absolute energies when false. If ordinary + ``IO.outputWF`` is requested without an explicit band list and both limits + remain zero, the default relative range selects no range and CONQUEST + disables wavefunction output with a warning. + + *default*: T + +IO.write_proj_DOS (*boolean*) + Writes the ordinary and overlap-scaled wavefunction coefficients required + for projected-DOS post-processing. This option also enables wavefunction + output and is available only for static diagonalisation calculations. + + *default*: F IO.TimingOn (*boolean*) Whether time information will be measured and written to output @@ -449,19 +504,19 @@ minE.SCTolerance (*real*) minE.SupportVariations (*integer*) Maximum number of support-function iterations - *default*: 20 + *default*: 10 -minE.PreconditionBlips(*boolean*) +minE.PreconditionBlips (*boolean*) Should blip variation be pre-conditioned? Pre-conditioning is (at present) more memory-intensive than it should be, but is efficient - *default*: F + *default*: T when ``minE.VaryBasis`` is true; otherwise disabled minE.GlobalTolerance (*boolean*) Are the convergence criteria applied to minimisation summed over the whole system, or per atom? - *default*: T + *default*: T for molecular dynamics; F for other run types Go to :ref:`top `. @@ -488,7 +543,7 @@ SC.LinearMixingFactor\_SpinDown (*real*) SC.LinearMixingEnd (*real*) Tolerance for end of Pulay mixing - *default*: self-consistency tolerance + *default*: ``minE.SCTolerance`` multiplied by :math:`10^{-4}` SC.LateStageReset (*integer*) If using GR-Pulay, how often is residual calculated fully (rather than interpolated) ? @@ -510,18 +565,6 @@ SC.MaxPulay (*integer*) *default*: 5 -SC.ReadAtomicDensityFile (*string*) - Filename for radial tables of atomic density (*rarely* used: normally generated from PAOs) - - default: - -SC.AtomicDensityFlag (*string*) - values: pao/read - - Flag determining how atomic densities should be found - - *default*: pao - SC.KerkerPreCondition (*boolean*) Flag determining if Kerker precondition is to be used. @@ -547,7 +590,9 @@ SC.MetricFactor (*real*) SC.MakeInitialChargeFromK (*boolean*) Flag determining whether initial charge is made from the density matrix - *default*: T + This option is disabled if ``General.LoadDM`` is false. + + *default*: T when ``General.LoadDM`` is true; F otherwise Go to :ref:`top `. @@ -566,6 +611,34 @@ DM.SolutionMethod (*string*) *default*: diagon +DM.DFTplusU (*boolean*) + Enables the DFT+U correction. A ``DFTplusU`` block is required when this + option is enabled. DFT+U is not available with a blip basis; CONQUEST + disables it with a warning in that case. + + *default*: F + +DM.WriteOccMat (*boolean*) + Prints the DFT+U occupation matrices in the main output. This option is + read only when ``DM.DFTplusU`` is true. + + *default*: T when ``DM.DFTplusU`` is true + +DFTplusU (*block*) + Defines the PAO projector and Hubbard U value for each corrected species. + Each line has the form: + + :: + + species n l zeta U + + ``species`` is the species number, ``n`` and ``l`` are the principal and + angular-momentum quantum numbers, ``zeta`` selects the zeta function, and + ``U`` is given in Hartree. The current implementation stores one projector + definition per species, so each corrected species should be listed once. + + *default*: none; required when ``DM.DFTplusU`` is true + DM.L\_range (*real*) Cutoff applied to L matrix (total energy will converge with increasing range; suggested minimum for O(N) calculations is twice largest support function range; @@ -628,11 +701,12 @@ DM.InvSDeltaOmegaTolerance (*real*) *default*: 0.0001 -DM.ConstantMu (*boolean*) - Switches between fixed Fermi level (T) and fixed number of electrons (F). You - are *strongly* recommended to leave at default +DM.ConstantMu (*presence flag*) + Fix the chemical potential at ``DM.mu``. Any occurrence of this tag, + including ``DM.ConstantMu F``, enables fixed-chemical-potential mode. Omit + the tag for the normal fixed-electron-number mode. - *default*: F + *default*: absent (fixed number of electrons) DM.mu (*real*) Value of Fermi level for fixed Fermi level calculations @@ -647,22 +721,43 @@ Diagonalisation --------------- Diag.NumKpts (*integer*) - Number of all k-points. No symmetry is applied. + Number of explicitly listed k-points, or the number of generated points on + each line when ``Diag.KspaceLines`` is true. No symmetry is applied. - *default*: + *default*: 1 for explicit points; 2 points per line -Diag.Kpoints (*block*)  - Lists fractional coordinates and weights of all k-points: ``x_fract y_fract z_fract weight`` - Generates the Monkhorst-Pack mesh, an equally spaced mesh of k-points. +Diag.Kpoints (*block*) + Lists fractional coordinates and weights of explicit k-points as + ``x_fract y_fract z_fract weight``. If this block is absent and neither a + Monkhorst-Pack mesh nor k-space lines are requested, CONQUEST uses the + Gamma point with unit weight. - *default*: + *default*: absent (Gamma point only) Diag.MPMesh (*boolean*) - Switches on/off the Monkhorst-Pack mesh. Note: if this keyword is present in - the input file, the keyword **Diag.NumKpts** and the block **Kpoints** will - be ignored. + Switches on/off the Monkhorst-Pack mesh. When true, explicit k-points and + k-space lines are ignored. + + *default*: F + +Diag.KspaceLines (*boolean*) + Generates k-points along lines in reciprocal space. This option is used + only when ``Diag.MPMesh`` is false. See the + :ref:`band-structure procedure ` for an example. - *default*: + *default*: F + +Diag.NumKptLines (*integer*) + Number of k-space lines used when ``Diag.KspaceLines`` is true. + + *default*: 1 + +Diag.KpointLines (*block*) + Lists the start and end points of each k-space line in fractional + coordinates. The block must contain two entries for every line specified + by ``Diag.NumKptLines``. + + *default*: none; required when ``Diag.KspaceLines`` is true Diag.MPMesh[X/Y/Z] (*integer*) Specifies the number n of k-points along the x(y,z) axis. @@ -692,11 +787,13 @@ Diag.PaddingHmatrix (*boolean*) Diag.BlockSizeR (*integer*) Block size for rows (See next). - From v1.4, the default value is 32 when Diag.PaddingHmatrix is true. + The default value is 32 when ``Diag.PaddingHmatrix`` is true. When padding + is disabled and this keyword is omitted, the value is determined + automatically. It is recommended to check the efficiency (CPU time) on your platform by changing this value. Usually 20-40 is appropriate. - *default*: 32 or Determined automatically (if Diag.PaddingHmatrix= true) + *default*: 32 with padding; otherwise determined automatically Diag.BlockSizeC (*integer*) R ... rows, C ... columns @@ -709,7 +806,12 @@ Diag.BlockSizeC (*integer*) If Diag.PaddingHmatrix is set to true then the block sizes can take any value, but BlockSizeR and BlockSizeC must be the same. - *default*: Determined automatically + With padding disabled, omitting ``Diag.BlockSizeR`` causes both block sizes + to be determined automatically. If ``Diag.BlockSizeR`` is specified, + ``Diag.BlockSizeC`` is read separately and defaults to 1 if omitted. + + *default*: equal to ``Diag.BlockSizeR`` with padding; otherwise as described + above Diag.MPShift[X/Y/Z] (*real*) Specifies the shift *s* of k-points along the x(y,z) axis, in fractional @@ -733,13 +835,16 @@ Diag.kT (*real*) *default*: 0.001 -Diag.IntegerOccs (*logical*) +Diag.IntegerOccs (*boolean*) Specifies integer occupation of eigenstates (overrides ``Diag.SmearingType``) *default*: False -Diag.AdjustEf (*logical*) - Turns off the automatic adjustment of Ef to be mid-gap (unlikely to be needed). +Diag.AdjustEf (*boolean*) + With Fermi-Dirac occupations, move the Fermi level to the middle of a + detected band gap. Set this false to disable the adjustment. This option + does not affect Methfessel-Paxton or integer occupations; integer + occupations always use the midgap Fermi level. *default*: True @@ -829,7 +934,7 @@ AtomMove.NumSteps (*integer*) AtomMove.MaxForceTol (*real*) The structure optimisation will stop when the maximum force component is less - than **MD.MaxForceTol** + than ``AtomMove.MaxForceTol``. *default*: 0.0005 Ha/bohr @@ -851,13 +956,10 @@ AtomMove.IonTemperature (*real*) *default*: 300 K for MD, 0 for Quench MD or FIRE AtomMove.ReadVelocity (*boolean*) - Read velocity from file ``md.checkpoint`` (when ``AtomMove.RestartRun T``) - - or ``velocity.dat`` (when ``AtomMove.RestartRun F``, very rare) + Read velocities from ``md.checkpoint`` during an MD restart, or from + ``velocity.dat`` otherwise - *default*: F (when ``AtomMove.RestartRun F``) - - or T (when ``AtomMove.RestartRun T``) + *default*: T when ``AtomMove.RestartRun`` is true; F otherwise AtomMove.AppendCoords (*boolean*) Chooses whether to append coordinates to ``UpdatedAtoms.dat`` during atomic @@ -866,31 +968,37 @@ AtomMove.AppendCoords (*boolean*) *default*: T AtomMove.OutputFreq (*integer*) - Frequency of output of information. *Not properly implemented* + Frequency, in ionic steps, for writing MD frames to ``md.frames`` and, when + ``MD.TDEP`` is true, TDEP-compatible lattice-dynamics files. It also + supplies the default MD frequency for the XSF and extended-XYZ trajectory + files. *default*: 50 AtomMove.WriteXSF (*boolean*) - Write atomic coordinates to ``trajectory.xsf`` for ``AtomMove.TypeOfRun = md`` or ``cg``, - every ``AtomMove.XsfFreq`` steps + Write atomic coordinates to ``trajectory.xsf``. During MD, frames are + written every ``AtomMove.XsfFreq`` steps. Structure-relaxation methods write + a frame after every ionic step. *default*: T AtomMove.XsfFreq (*integer*) - Frequency of output of atomic coordinates to ``trajectory.xsf`` + Frequency, in MD steps, of atomic-coordinate output to ``trajectory.xsf``. + This setting is not used by structure-relaxation methods. *default*: same as ``AtomMove.OutputFreq`` -AtomMove.WriteXYZ (*boolean*) - Write atomic coordinates to ``trajectory.xyz`` for ``AtomMove.TypeOfRun = md``, - every ``AtomMove.XyzFreq`` steps +AtomMove.WriteExtXYZ (*boolean*) + Write coordinates, lattice vectors, energy, forces and stress in extended + XYZ format to ``trajectory.xyz``. Moving-atom runs write frames every + ``AtomMove.XyzFreq`` steps; a static run writes one frame. - *default*: T + *default*: F AtomMove.XyzFreq (*integer*) - Frequency of output of atomic coordinates to ``trajectory.xyz`` + Frequency of output to ``trajectory.xyz`` during moving-atom runs. - *default*: same as ``AtomMove.OutputFreq`` + *default*: same as ``AtomMove.OutputFreq`` for MD; 1 for other run types AtomMove.TestForces (*boolean*) Flag for testing forces with comparison of analytic and numerical calculations. @@ -901,10 +1009,10 @@ AtomMove.TestForces (*boolean*) AtomMove.TestAllForces (*boolean*) Switch to test *all* force contributions or not - *default*: F + *default*: T AtomMove.CalcStress (*boolean*) - Toggle calculation of the stress tensor. Switching off can improve performace. + Toggle calculation of the stress tensor. Switching off can improve performance. *default*: T @@ -983,6 +1091,8 @@ AtomMove.OptCell.Constraint (*string*) ``volume``: minimize the total energy by scaling each simulation cell dimension by the same global scaling factor. Search directions are set by the mean stress. + *default*: none + AtomMove.TestSpecificForce (*integer*) Label for which force contribution to test. Note that for PAOs non-local Pulay and Hellman-Feynman forces are found together as part of the HF calculation; @@ -1016,11 +1126,12 @@ AtomMove.TestForceDelta (*real*) *default*: 10\ :math:`^{-5}` bohr AtomMove.RestartRun (*boolean*) - Restart a MD run. Note that this will set ``General.LoadL T``, - ``AtomMove.MakeInitialChargeFromSC T`` and ``XL.LoadX T`` if using the - extended Lagrangian. The atomic coordinates will be read from - ``md.positions`` and the velocities and extended system variables from - ``md.checkpoint``. + Restart an MD run. Unless explicitly overridden, this defaults + ``General.LoadDM``, ``SC.MakeInitialChargeFromK`` and + ``AtomMove.ReadVelocity`` to T. It also defaults ``XL.LoadX`` to T for + XL-BOMD, and ``Basis.LoadCoeffs`` to T for MSSF and blip-basis calculations. + Atomic coordinates are read from ``md.position``; velocities and extended- + system variables are read from ``md.checkpoint``. *default*: F @@ -1049,7 +1160,7 @@ AtomMove.SkipEarlyDM (*boolean*) AtomMove.McWeenyFreq (*integer*) McWeeny step is applied every N steps (with “AtomMove.ReuseDM T”) - *default*: + *default*: 0 (disabled) AtomMove.ExtendedLagrangian (*boolean*) Selects XL-BOMD (with “AtomMove.ReuseDM T”) @@ -1059,7 +1170,8 @@ AtomMove.ExtendedLagrangian (*boolean*) AtomMove.FixCentreOfMass (*boolean*) Remove the centre of mass velocity at every time step - *default*: T + *default*: T for MD, except F for FIRE; disabled when any atomic coordinate + is fixed and unused for other run types Go to :ref:`top `. @@ -1069,52 +1181,71 @@ Molecular Dynamics ------------------ MD.Ensemble (*string*) - values: nve/nvt/npt/nph + values: nve/nvt/npt The molecular dynamics ensemble *default*: nve MD.Thermostat (*string*) - values: none/nhc/berendsen/svr + values: none/nhc/svr Thermostat type ``none`` No thermostat (used for calculating temperature only) - ``berendsen`` - Berendsen weak coupling thermostat + ``nhc`` + Nose-Hoover chain thermostat ``svr`` Stochastic velocity rescaling - *default*: none + *default*: none for NVE; nhc for NVT and NPT MD.Barostat (*string*) - values: none/berendsen/iso-mttk/ortho-mttk/mttk + values: none/pr/mttk - Barostat type. The following are the only valid thermostat/barostat - combinations for the NPT ensemble: ``berendsen``/ ``berendsen``, - ``nhc``/ ``pr``, ``svr``/ ``pr`` + Barostat type. ``pr`` selects the Parrinello-Rahman implementation; + ``mttk`` selects the Martyna-Tobias-Tuckerman-Klein implementation. ``none`` No barostat (used for calculating pressure only) - ``berendsen`` - Berendsen weak coupling barostat ``pr`` Parrinello-Rahman (extended system) barostat + ``mttk`` + Martyna-Tobias-Tuckerman-Klein extended-system barostat - *default*: none + *default*: none for NVE and NVT; pr for NPT + +MD.CellConstraint (*string*) + values: fixed/volume/xyz + + Controls the cell degrees of freedom used by the barostat. ``volume`` uses + one isotropic cell degree of freedom, preserving the cell shape. ``xyz`` + allows the three orthorhombic cell lengths to change independently and is + implemented for the ``pr`` barostat. The ``mttk`` implementation is + isotropic and should be used with ``volume``. When no barostat is active, + CONQUEST sets the constraint to ``fixed``; use ``MD.Barostat none`` for + fixed-cell dynamics. + + *default*: volume; changed to fixed when no barostat is active + +MD.EquilSteps (*integer*) + Number of initial equilibration steps. A positive value is retained only + when ``MD.Thermostat`` is ``svr``; otherwise CONQUEST issues a warning and + resets it to zero. + + *default*: 0 MD.tauT (*real*) - Coupling time constant for thermostat. Required for Berendsen thermostat, or - if ``MD.CalculateXLMass = T``. Note that this number means different things - for the Berendsen and NHC thermostats. + Thermostat coupling time in fs. For SVR this is a relaxation timescale; for + NHC it sets the thermostat frequency when ``MD.CalculateXLMass`` is true. - *default*: 1.0 + *default*: 50 fs for SVR; 10 times ``AtomMove.Timestep`` for NHC MD.TDrag (*real*) - Add a drag coefficient to the thermostat. The thermostat velocities are - reduced by a factor :math:`1 - \tau/D_T` every step. + Dimensionless strength of ad hoc damping applied to the ionic + Nose-Hoover-chain velocities. The internal damping factor also depends on + the timestep, ``MD.tauT``, ``MD.nMTS`` and ``MD.nYoshida``. *default*: 0.0 @@ -1128,7 +1259,7 @@ MD.CellNHC (*boolean*) *default*: T -MD.NHCMass (*blocks*) +MD.NHCMass (*block*) :math:` \ldots` Masses of NHC heat baths @@ -1141,22 +1272,22 @@ MD.CellNHCMass (*block*) *default*: 1 1 1 1 1 MD.BulkModulusEst (*real*) - Bulk modulus estimate for system. Only necessary for Berendsen weak pressure - coupling (``MD.Barostat = berendsen`` or ``MD.BerendsenEquil > 0``) + Reserved bulk-modulus estimate. Version 1.6 reads and stores this value, but + the supported barostat algorithms do not use it. *default*: 100 MD.tauP (*real*) - Coupling time constant for barostat. Required for Berendsen barostat, or if - MD.CalculateXLMass = T. Note that this number means different things for the - Berendsen and Parrinello-Rahman barostats. + Barostat coupling time in fs. It sets the barostat frequency when + ``MD.CalculateXLMass`` is true. - *default*: 10.0 (Berendsen) or 100.0 (MTTK) + *default*: 100 times ``AtomMove.Timestep`` MD.PDrag (*real*) - Add a drag coefficient to the barostat. The barostat velocities are - reduced by a factor :math:`1 - \tau/D_P` every step. This is useful - when the lattice parameters are varying rapidly. + Dimensionless strength of ad hoc damping applied to the cell/barostat + velocities and, when ``MD.CellNHC`` is true, the separate cell + Nose-Hoover-chain velocities. The internal damping factor also depends on + the timestep, ``MD.tauP``, ``MD.nMTS`` and ``MD.nYoshida``. *default*: 0.0 @@ -1183,11 +1314,6 @@ MD.nMTS (*integer*) *default*: 1 -MD.BerendsenEquil (*integer*) - Equilibrate the system for :math:`n` steps using Berendsen weak coupling - - *default*: 0 - MD.TDEP (*boolean*) Dump data in a format readable by the Temperature Dependent Effective Potential (TDEP) code. @@ -1220,12 +1346,12 @@ MD.VariableTemperatureRate (*real*) *default*: 0.0 -MD.InitialTemperature(*real*) +MD.InitialTemperature (*real*) Initial temperature. *default*: same as AtomMove.IonTemperature -MD.FinalTemperature(*real*) +MD.FinalTemperature (*real*) Final temperature. *default*: same as AtomMove.IonTemperature @@ -1243,7 +1369,8 @@ Spin.SpinPolarised (*boolean*) *default*: F Spin.FixSpin (*boolean*) - Determines if spin populations are to be fixed. Only read if **Spin.FixPolarised** is set. + Determines if spin populations are to be fixed. Only used if + ``Spin.SpinPolarised`` is true. *default*: F @@ -1264,65 +1391,71 @@ Go to :ref:`top `. DeltaSCF -------- -flag\_DeltaSCF (*boolean*) +minE.DeltaSCF (*boolean*) Selects delta SCF calculation - *default*: + *default*: F DeltaSCF.SourceLevel (*integer*) Eigenstate number to remove electron from (source) - *default*: + *default*: 0 DeltaSCF.TargetLevel (*integer*) Eigenstate number to promote electron to (target) - *default*: + *default*: 0 DeltaSCF.SourceChannel (*integer*) Spin channel for electron source - *default*: + *default*: 1 DeltaSCF.TargetChannel (*integer*) Spin channel for electron target - *default*: + *default*: 1 DeltaSCF.SourceNFold (*integer*) Allows selection of more than one level for excitation source (N-fold) - *default*: + *default*: 1 DeltaSCF.TargetNFold (*integer*) Multiplicity of target (N-fold) - *default*: + *default*: 1 DeltaSCF.LocalExcitation (*boolean*) - Select an excitation localised on a group of atoms + Select an excitation localised on a group of atoms. When enabled, the + ``cDFT.AtomGroups`` block and every atom-group block named within it are + required; CONQUEST stops if any of these blocks is absent. At least one of + ``DeltaSCF.HOMOLimit`` and ``DeltaSCF.LUMOLimit`` must be nonzero for the + localised-excitation analysis to remain enabled. - *default*: + *default*: F DeltaSCF.HOMOLimit (*integer*) How many states down from HOMO to search for localised excitation - *default*: + *default*: 0 DeltaSCF.LUMOLimit (*integer*) How many states up from LUMO to search for localised excitation - *default*: + *default*: 0 DeltaSCF.HOMOThresh (*real*) - (*please fill in*) + Threshold for identifying a localised occupied source state (sum of the + squared moduli of coefficients on the selected atoms) - *default*: + *default*: 0.5 DeltaSCF.LUMOThresh (*real*) - Threshold for identifying localised excitation (sum over square moduli of coefficients) + Threshold for identifying a localised unoccupied target state (sum of the + squared moduli of coefficients on the selected atoms) - *default*: + *default*: 0.5 Go to :ref:`top `. @@ -1334,34 +1467,34 @@ Constrained DFT (cDFT) cDFT.Perform\_cDFT (*boolean*) Selects cDFT operation - *default*: + *default*: F cDFT.Type (*integer*) values: 1 or 2 Selects constraint to be for absolute charge on groups (1) or difference between two groups (2) - *default*: + *default*: 2 cDFT.MaxIterations (*integer*) Maximum iterations permitted - *default*: + *default*: 50 cDFT.Tolerance (*real*) Tolerance on charge - *default*: + *default*: 0.001 cDFT.NumberAtomGroups (*integer*) Number of groups of atoms - *default*: + *default*: 1 cDFT.AtomGroups (*block*) - Block with each line specifying: Number of atoms, target charge, label for - block. For each line, there should be a corresponding block with the appropriate - label; the block consists of a list of atom numbers for the atoms in the group + Each line specifies a group index, number of atoms, target charge and block + label. For each line, a block with that label must list the atom numbers in + the group. Go to :ref:`top `. @@ -1406,10 +1539,10 @@ Go to :ref:`top `. vdW-DF ------ -vdWDFT.LDAFunctionalType (*string*) +vdWDFT.LDAFunctionalType (*integer*) Selects LDA functional to use with vdW-DF - *default*: + *default*: 3 Go to :ref:`top `. @@ -1421,7 +1554,7 @@ DFT-D2 DFT-D2\_range (*real*) DFT-D2 cutoff range (bohr) - *default*: + *default*: 23.0 Go to :ref:`top `. @@ -1448,7 +1581,7 @@ XL.PropagateL (*boolean*) XL.Dissipation (*boolean*) Selects the addition of dissipative force - *default*: + *default*: F XL.MaxDissipation (*integer*) Order of dissipative force term @@ -1474,8 +1607,8 @@ Advanced and obscure tags .. _advanced_general_tags: -General -******* +General (advanced) +****************** General.LoadInvS (*boolean*) Selects loading of inverse S matrix from previous step (not @@ -1507,11 +1640,6 @@ General.EwaldAccuracy (*real*) *default*:1\ :math:`\times`\ 10\ :math:`^{-10}` -General.CheckDFT (*boolean*) - Calculates DFT energy using output density - - *default*: F - General.AverageAtomicDiameter (*real*) Related to space-filling @@ -1525,6 +1653,8 @@ General.GapThreshold (*real*) General.only_Dispersion (*boolean*) Selects only DFT\_D2 calculation (no electronic structure etc) + *default*: F + General.MixXCGGAInOut (*real*) For non-SCF calculations only, chooses how to mix the proportions of GGA XC stress contribution (from the change of the electron density @@ -1539,8 +1669,8 @@ Go to :ref:`top `. .. _advanced_atomic_spec_tags: -Atomic Specification -******************** +Atomic Specification (advanced) +******************************* Atom.ValenceCharge (*real*) Valence charge of species (e.g. 4 for carbon, 6 for oxygen) @@ -1652,7 +1782,7 @@ IO.PdbAltLoc (*string*) IO.PdbOut (*boolean*) Format of the output coordinate file. Writes a PDB file if set to T. In that case, either the input must be in pdb format or a PDB “template” file needs to - be specified (keyword General.PdbTemplate) + be specified with ``IO.PdbTemplate``. *default*: F @@ -1663,7 +1793,7 @@ IO.PdbTemplate (*string*) it will also be used as the template, although this can still be overwritten with this keyword - *default*: coordinate file + *default*: ``IO.Coordinates`` when ``IO.PdbIn`` is true; no template otherwise IO.AtomOutputThreshold (*integer*) Threshold below which atomic positions are output on @@ -1688,23 +1818,16 @@ Basis.BasisSet (*string*) - PAOs — Pseudo-atomic orbitals :cite:`e-Artacho1999` - - blips (default) — B-splines :cite:`e-Hernandez1997` + - blips — B-splines :cite:`e-Hernandez1997` *default*: PAOs -Basis.LoadBlip (*boolean*) - Load blip or PAO coefficients from file. If set to T, for blips the code will - look for a set of files containing blip coefficients, which is taken to be - ``blip_coeffs.nnn``, where ``nnn`` is processor number (padded with zeroes); - for PAOs, the code will look for a *single* file which is ``supp_pao.dat`` - by default, but can be set with ``Basis.SupportPaoFile`` - - *default*: F - -Basis.SupportPaoFile (*string*) - Specifies filename for PAO coefficients +Basis.LoadCoeffs (*boolean*) + Load existing support-function coefficients. Blip calculations read the + per-process ``blip_coeffs`` files; MSSF calculations read the dumped + ``SFcoeff`` matrix files. - *default*: ``supp_pao.dat`` + *default*: F, or T when ``AtomMove.RestartRun`` is true (for blips and MSSF) Basis.UsePulayForPAOs (*boolean*) Determines whether to use Pulay DIIS for minimisation of PAO basis coefficients @@ -1770,17 +1893,12 @@ Basis.TestBasisGrad\_H (*boolean*) *default*: F -Basis.PAOs\_OneToOne (*boolean*) - Assign PAOs to individual support functions (implies no support function optimisation) - - *default*: F - Go to :ref:`top `. .. _advanced_grid_tags: -Integration Grid -**************** +Integration Grid (advanced) +*************************** Grid.PointsAlong[X/Y/Z] (*integer*) Grid points along x (y,z). Overwrites the values set by **Grid.GridCutoff**. @@ -1810,4 +1928,3 @@ Go to :ref:`top `. :style: unsrt Go to :ref:`top `. - diff --git a/docs/installing.rst b/docs/installing.rst index 621f148c2..f0ae6a908 100644 --- a/docs/installing.rst +++ b/docs/installing.rst @@ -88,14 +88,14 @@ Compiler flags to enable OpenMP are dependent on the vendor, but should be speci On some systems, the default stack size for OpenMP is set to be rather small, and this can cause a segmentation fault when running with multiple threads. We recommend testing the effect of the environment variable ``OMP_STACKSIZE`` (and suggest setting it to 50M or larger as a first test). -Most OpenMP multi-threading in CONQUEST uses the ``runtime`` schedule. This means the type of scheduling of work to threads can be set by the user by setting the ``OMP_SCHEDULE`` `variable`_. If the variable is unset, OpenMP will use a default implementation defined schedule. +Most OpenMP multi-threading in CONQUEST uses the ``runtime`` schedule. This means the type of scheduling of work to threads can be set by the user by setting the ``OMP_SCHEDULE`` `variable `_. If the variable is unset, OpenMP will use a default implementation defined schedule. Go to :ref:`top ` .. _install_spack: Installing with Spack ------------ +--------------------- CONQUEST and all of its dependencies can be installed with `Spack `_. The CONQUEST package requires Spack v0.21 or later. If Spack isn't available or up to date on your @@ -117,14 +117,14 @@ The build can be customized by adding options to the The CONQUEST package includes variants for OpenMP support and different matrix multiplication kernels; more details can be found in the `Spack CONQUEST package `_. Installing on Ubuntu ------------ +-------------------- CONQUEST can be compiled on Ubuntu after installing the required packages. The below instructions are given for Ubuntu 22.04 LTS and Ubuntu 24.04 LTS. The source files will be downloaded into the ``${USER}/local/src`` directory. The ${USER} variable will be automatically replaced by the current username. If compilation is successful, the executable file can be found in ``${USER}/local/src/conquest_master/bin``. Install needed packages -~~~~~~~~~~~~~~~~ +~~~~~~~~~~~~~~~~~~~~~~~ .. code-block:: bash @@ -218,7 +218,7 @@ Prepare makefile EOF Compile CONQUEST -~~~~~~~~~~~~~~~ +~~~~~~~~~~~~~~~~ .. code-block:: bash diff --git a/docs/moldyn.rst b/docs/moldyn.rst index b57ef7940..d26331930 100644 --- a/docs/moldyn.rst +++ b/docs/moldyn.rst @@ -26,7 +26,7 @@ conserved quantity of the dynamics. Although the molecular dynamics integrators used in CONQUEST are time reversible, *the SCF procedure is not*. Therefore tight convergence (``minE.SCTolerance`` for diagonalisation, ``minE.LTolerance`` for linear scaling) is -necessary. In the case of diagonalisation, SCF tolerance of ``1E-6`` is +necessary. In the case of diagonalisation, an SCF tolerance of ``1E-6`` is typically enough to negate the drift. However, extended-Lagrangian Born-Oppenheimer MD (XL-BOMD) :cite:`md-Niklasson2008`, currently only implemented for O(N), essentially makes the SCF component of the MD @@ -34,7 +34,7 @@ time-reversible by adding the electronic degrees of freedom to the Lagrangian, relaxing the constraint on ``minE.LTolerance`` --- although it is still somewhat dependent on the ensemble. In the NVE and NVT ensembles, a L-tolerance of ``1E-5`` has been found to be -sufficient to give good energy conservations, decreasing to ``1E-6`` +sufficient to give good energy conservation, decreasing to ``1E-6`` in the NPT ensemble. The following flags are required for XL-BOMD: :: @@ -58,30 +58,36 @@ setting, This will do several things: it will read the atomic coordinates from ``md.position`` and read the ``md.checkpoint`` file, which contains the -velocities and extended system (Nose-Hoover chain and cell) variables. Depending +velocities and extended-system (Nose-Hoover chain and cell) variables. Depending on the value of ``DM.SolutionMethod``, it will read the K-matrix files -(``diagon``) or the L-matrix files (``ordern``), and if XL-BOMD is being used, -the X-matrix files. Finally, it will *append* new data to the ``md.stats`` and -``md.frames`` files, but it will overwrite all other files, including -``Conquest_out``. Note that this flag is equivalent to setting the following: +(``diagon``) or the L-matrix files (``ordern``). Finally, it will *append* new +data to the ``md.stats`` and ``md.frames`` files, but it will overwrite all +other files, including ``Conquest_out``. + +Unless explicitly overridden, ``AtomMove.RestartRun T`` changes the defaults of +the following keywords to ``T``: :: - General.LoadL T + General.LoadDM T SC.MakeInitialChargeFromK T - XL.LoadL T + AtomMove.ReadVelocity T + +For XL-BOMD it also defaults ``XL.LoadX`` to ``T``. For MSSF and blip-basis +calculations it defaults ``Basis.LoadCoeffs`` to ``T``, so that the saved +support-function coefficients are loaded. In addition to the files mentioned above, CONQUEST will try to read the K-matrix from ``Kmatrix2.i00.*`` when using diagonalisation or the L-matrix from ``Lmatrix2.i00.*`` when using O(N), and ``Xmatrix2.i0*.*`` if the extended-Lagrangian formalism is used. Note that metadata for these files is -stored in ``InfoGlobal.i00.dat`` which is also required when restarting. If the +stored in ``InfoGlobal.i00`` which is also required when restarting. If the calculation ended by hitting the walltime limit, the writing of these matrix -files may have been interrupted, rendering them unusable. In this case, the -calculation can be restarted by setting the above flags to ``F`` *after* setting -``AtomMove.RestartRun T``. Setting the flag ``General.MaxTime`` to some number -of seconds less (say 30 minutes) than the calculation wall time limit will force -the calculation to stop gracefully, preventing the aforementioned situation. +files may have been interrupted, rendering them unusable. In this case, set +``General.LoadDM F`` and, for XL-BOMD, ``XL.LoadX F`` *after* setting +``AtomMove.RestartRun T``. Setting ``General.MaxTime`` to some number of seconds +less (say 30 minutes) than the calculation wall time limit will force the +calculation to stop gracefully, preventing the aforementioned situation. Go to :ref:`top `. @@ -91,8 +97,9 @@ Visualising the trajectory -------------------------- Setting the flag ``AtomMove.WriteXSF T`` dumps the coordinates to the file -``trajectory.xsf`` every ``AtomMove.OutputFreq`` steps. The .xsf file can be -read using `VMD `_. A small VMD script, +``trajectory.xsf`` every ``AtomMove.XsfFreq`` steps (which defaults to +``AtomMove.OutputFreq``). The .xsf file can be read using +`VMD `_. A small VMD script, ``view.vmd`` is included with the code, and can be invoked using, ``vmd -e view.vmd`` @@ -106,7 +113,7 @@ Go to :ref:`top `. TDEP output ----------- -CONQUEST molecular dynamics data can be used to perform lattice dyanmical +CONQUEST molecular dynamics data can be used to perform lattice dynamical calculations using the `Temperature Dependent Effective Potential (TDEP) `_ code. Setting the flag ``MD.TDEP T`` will make conquest dump configurations, forces and metadata in a format @@ -187,11 +194,9 @@ Go to :ref:`top `. Isobaric-Isothermal (NPT) ensemble ++++++++++++++++++++++++++++++++++ -There is one implemented barostat at present, the extended -system, Parrinello-Rahman :cite:`md-Parrinello1981`. At present the -barostat should be treated as a beta-version implementation, which -will be fully characterised and made robust for the full release of -the code. +Two extended-system barostats are implemented: Parrinello-Rahman +:cite:`md-Parrinello1981` and Martyna-Tobias-Tuckerman-Klein (MTTK) +:cite:`t-Martyna1996`. 1. Parrinello-Rahman @@ -218,20 +223,47 @@ NHC tend to be more severe due to coupling of the cell and atomic motions. They are dependent on the system, so it is advised that you find a combination of these parameters that gives the best energy conservation. The cell is thermostatted using a separate Nose-Hoover chain to the atoms by default, but -they can be controlled with the same chain by setting ``MD.CellNHC F``. An *ad -hoc* drag factor specified by ``MD.PDrag`` reduces the thermostat and cell -velocities at every timestep to damp out the ringing fluctuations. In this case, -they are reduced by :math:`10/200 \simeq 5\%`, which strictly speaking breaks the NPT -dynamics, but not significantly, and the stability is significantly improved. +they can be controlled with the same chain by setting ``MD.CellNHC F``. The +*ad hoc* damping controlled by ``MD.PDrag`` acts on the cell/barostat velocities +and, when present, the separate cell-thermostat velocities. Each application +uses the factor + +.. math:: + + 1 - \frac{D_P\,\Delta t}{\tau_P n_{\mathrm{MTS}} n_{\mathrm{YS}}}, + +where :math:`D_P` is ``MD.PDrag``. This factor is applied within the integration +substeps, potentially more than once per full MD timestep. A nonzero drag +perturbs the formal NPT dynamics, but can help damp ringing fluctuations. Note that the NPT ensemble can also be generated correctly by thermostatting using the SVR thermostat, although the meaning of the parameter ``MD.tauT`` is different in this case, as in NVT dynamics. +2. Martyna-Tobias-Tuckerman-Klein + +The MTTK barostat is implemented for isotropic volume fluctuations and is +coupled to the Nose-Hoover-chain thermostat: + +:: + + AtomMove.IonTemperature 300.0 + AtomMove.TargetPressure 10.0 + MD.Ensemble npt + MD.Thermostat nhc + MD.Barostat mttk + MD.CellConstraint volume + MD.tauT 100 + MD.tauP 200 + +Use the Parrinello-Rahman barostat instead when the three orthorhombic cell +lengths must vary independently with ``MD.CellConstraint xyz``. + Postprocessing tools -------------------- -Details of Python post-processing tools for CONQUEST can be found in :ref:`et_md_scripts`. +The Python utility supplied for analysing CONQUEST molecular-dynamics output is +described in :ref:`et_md_scripts`. Go to :ref:`top `. diff --git a/docs/post-proc.rst b/docs/post-proc.rst index e3c6b0f34..2d72e71c8 100644 --- a/docs/post-proc.rst +++ b/docs/post-proc.rst @@ -31,9 +31,9 @@ produced during the CONQUEST run are available for the post-processing flags from the CONQUEST run that generated the output, and some utility-specific flags that are detailed below. -**Note also** that projected DOS, band density and STM simulation are -not at present compatible with multi-site support functions (MSSF), -though we hope to implement this soon. +Projected DOS, band density and STM simulation are available for calculations +using primitive PAOs or multi-site support functions (MSSFs), but not for +calculations using blip functions. Go to :ref:`top `. @@ -185,8 +185,8 @@ Density of states (DOS) +++++++++++++++++++++++ Setting ``Process.Job dos`` will produce a total density of states -(DOS) for the system, using the eigenvalues output by CONQUEST. The -following parameters can be set: +(DOS) for the system, using the ``eigenvalues.dat`` file output by a +CONQUEST diagonalisation calculation. The following parameters can be set: :: @@ -238,21 +238,17 @@ given state, :math:`n`, onto an atom :math:`i` can be written as S_{i\alpha,j\beta}c^{n\mathbf{k}}_{j\beta}`. The projected DOS is constructed using these projections. -If using :ref:`pseudo-atomic orbitals (PAOs) ` as the -basis set, then the atom-projected DOS can be further resolved by -angular momentum (either just :math:`l` or both :math:`l` and -:math:`m`). If using :ref:`pseudo-atomic orbitals (PAOs) -` with :ref:`multi-site support functions ` or -:ref:`blip functions ` then it is not possible to -decompose the DOS any further (in future, it may be possible to -resolve the MSSF coefficients into the individual PAOs, and hence -decompose pDOS by angular momentum). To output the necessary -coefficients to produce atom-projected DOS, a CONQUEST run must be -performed with the following parameters set: +For calculations using :ref:`pseudo-atomic orbitals (PAOs) `, +either directly or through :ref:`multi-site support functions `, +the atom-projected DOS can be further resolved by angular momentum (either +just :math:`l` or both :math:`l` and :math:`m`). For MSSF calculations, +CONQUEST outputs the required coefficients in the underlying PAO basis. +Projected-DOS post-processing is not available for calculations using +:ref:`blip functions `. To output the necessary coefficients, +a static diagonalisation calculation must be performed with: :: - IO.writeDOS T IO.write_proj_DOS T As for the DOS, very high Brillouin zone sampling is required for @@ -287,10 +283,9 @@ following flags can be set: Process.pDOS_lm_resolved T Note that only one of these is needed, depending on what level of -resolution is required. At present, angular momentum resolution is -only available for the PAO basis set (not MSSF or blips) though it -is under development for the MSSF basis (by projection onto the -underlying PAO basis). +resolution is required. Angular-momentum resolution is available for both +primitive-PAO and MSSF calculations, but not for calculations using blip +functions. The energy range for the projected DOS can also be specified: @@ -322,9 +317,9 @@ Go to :ref:`top `. Rotated pDOS ------------ -CONQUEST, by default, projects the orbitals along the simulation cell axes. Sometimes it is convenient to be able to project onto a coordinate system defined by an atom's local environment, e.g. bonds, to facilitate chemical analysis. CONQUEST supports rotating the wavefunction coefficients either by inputting the basis of the final coordinate system, using the algorithm implemented by Maintz *et al* and Romanowski *et al* [:cite:`pp-maintz2016`, :cite:`pp-romanowski2008`] or active Euler angles in the extrinsic :math:`zyz` convention. +CONQUEST, by default, projects the orbitals along the simulation cell axes. Sometimes it is convenient to be able to project onto a coordinate system defined by an atom's local environment, e.g. bonds, to facilitate chemical analysis. CONQUEST supports rotating the wavefunction coefficients either by inputting the basis of the final coordinate system, using the algorithms described by Maintz *et al.* and Romanowski *et al.* :cite:`pp-maintz2016,pp-romanowski2008`, or by supplying active Euler angles in the extrinsic :math:`zyz` convention. -There are 4 operating modes that are possible: +There are four operating modes: #. Supply a set of axes to apply to a set of (or all) atoms (mode 0) #. Supply 3 Euler angles to apply to a set of (or all) atoms (mode 1) @@ -365,7 +360,15 @@ If Euler angles are desired, set ``Process.RotatePDOSMode 1`` and in ``Conquest_ %endblock pDOSEuler -Here, ``A``: the atom number to assign the Euler angles to. The default angle units is in degrees: ``Process.RotatePDOSAngle deg`` and units cannot be mixed. In the extrinsic :math:`zyz` convention, 3 consecutive rotations happen about a fixed set of coordinates, i.e. the CONQUEST cell simulation axes (as CONQUEST only supports orthorhombic cells), first about the :math:`z`-axis by :math:`\gamma` , then :math:`y`-axis by :math:`\beta` and then about :math:`z`-axis by :math:`\alpha` in a right-handed sense (i.e., looking down any of the axes towards the origin will mean the positive direction of rotation is anticlockwise). +Here, ``A`` is the atom number to assign the Euler angles to. The default angle +unit is degrees (``Process.RotatePDOSAngle deg``), and units cannot be mixed. In +the extrinsic :math:`zyz` convention, three consecutive rotations happen about +a fixed set of coordinates, i.e. the CONQUEST simulation-cell axes (CONQUEST +only supports orthorhombic cells): first about the :math:`z` axis by +:math:`\gamma`, then the :math:`y` axis by :math:`\beta`, and then the +:math:`z` axis by :math:`\alpha`, in a right-handed sense (looking down any of +the axes towards the origin, the positive direction of rotation is +anticlockwise). Determining local axes is possible depending on local geometries. @@ -381,9 +384,9 @@ Determining local axes is possible depending on local geometries. where ``A``, ``B``, ``C``, ``D`` are integers with different conditions: * ``A``: the atom number to perform the rotation for -* ``B``: if it is ``0``, assume square planar geometry, if ``1`` assume octahedral. This controls nearest-neighbour searching. This wil search for the nearest 4 or 6 neighbours respectively. +* ``B``: if it is ``0``, assume square planar geometry; if ``1``, assume octahedral. This controls nearest-neighbour searching for the nearest four or six neighbours, respectively. * ``C``: if ``-1``, set the principal direction along the shortest bond. If ``0``, set along the longest bond. If it is a positive integer, it is interpreted as an atom number which must be a neighbour of the atom specified in ``A``. All runs will output the atom neighbours, so setting this to ``0`` or ``-1`` as a first run is recommended unless the neighbours are known beforehand -* ``D``: if ``0``, the second direction is chosen by the bond which changes the least under projection to the plane defined by the bond vector calculated from the neighbour specified in ``C``. If ``D`` is a positive integer, then it is interpreted as an atom number which must be a neighbour of the atom specified in ``A``. The code will error out if it is the same as ``C``. A warning will be displayed if the direction towards neighbour ``D`` is determined to deviate significantly away from orthogonal. +* ``D``: if ``0``, the second direction is chosen by the bond which changes the least under projection to the plane defined by the bond vector calculated from the neighbour specified in ``C``. If ``D`` is a positive integer, then it is interpreted as an atom number which must be a neighbour of the atom specified in ``A``. A warning will be displayed if the direction towards neighbour ``D`` is determined to deviate significantly away from orthogonal. If the rotation axis and angle are known beforehand, using mode 3 can be convenient. @@ -405,7 +408,7 @@ For users interested in verbose output, set ``Process.RotatePDOSDebug T`` (defau Go to :ref:`top `. -.. _pp_band_str +.. _pp_band_str: Band structure ++++++++++++++ @@ -414,8 +417,8 @@ The band structure of a material can be generated by CONQUEST by performing a non-self-consistent calculation, after reading a well-converged charge density: set ``minE.SelfConsistent F`` and ``General.LoadRho T`` (remember that to write a converged charge density from CONQUEST you set ``IO.DumpChargeDensity T``). -The k-points required can be specified as lines of points in k-space; -setting ``Diag.KspaceLines T`` enables this (replacing the usual MP mesh), while the number of lines +The k-points required can be specified as lines of points in k-space. With +``Diag.MPMesh F``, setting ``Diag.KspaceLines T`` enables this, while the number of lines (e.g. Gamma to L; L to X; would be two lines) is set with ``Diag.NumKptLines`` and the number of points along a line with ``Diag.NumKpts``. The k-point lines themselves are set with a block labelled ``Diag.KpointLines`` which should have @@ -450,3 +453,8 @@ with ``Process.min_DOS_E`` and ``Process.max_DOS_E``. Go to :ref:`top `. +.. bibliography:: references.bib + :cited: + :labelprefix: PP + :keyprefix: pp- + :style: unsrt diff --git a/docs/strucrelax.rst b/docs/strucrelax.rst index 3b89d843e..32aa475f3 100644 --- a/docs/strucrelax.rst +++ b/docs/strucrelax.rst @@ -14,7 +14,7 @@ conjugate gradients at present, though L-BFGS will be implemented. Setting ``AtomMove.WriteXSF T`` for all flavours of optimisation will dump the trajectory to the file ``trajectory.xsf``, which can be visualised using `VMD -`_ and `XCrysDen `_. +`_ and `XCrySDen `_. Setting ``AtomMove.AppendCoords T`` will append the structure at each step to ``UpdatedAtoms.dat`` in the format of a CONQUEST structure input. @@ -82,12 +82,12 @@ For large initial forces or problematic cases where the relaxation algorithms fa downhill search direction, it may be worth trying quenched molecular dynamics, which propagates the equations of motion following a simple NVE approach, but resets the velocities to zero when the dot product of -force and velocity is zero. +force and velocity becomes negative. :: AtomMove.TypeOfRun md - AtomMove.QuenchedMD T + AtomMove.QuenchMD T AtomMove.MaxForceTol 5e-4 AtomMove.ReuseDM T diff --git a/docs/theory-md.rst b/docs/theory-md.rst index f818214e9..b0d9c7ee6 100644 --- a/docs/theory-md.rst +++ b/docs/theory-md.rst @@ -73,12 +73,12 @@ Go to :ref:`top `. Extended Lagrangian Born-Oppenheimer MD (XL-BOMD) ------------------------------------------------- -If the electronic density from the previous ionic step is used as an initila +If the electronic density from the previous ionic step is used as an initial guess for the next SCF cycle, a problem arises because this process breaks the time-reversibility of the dynamics. This is manifested as a gradual drift in the total energy in the case of a NVE simulation, or the conserved quantity in the case of non-Hamiltonian dynamics. The solution proposed by Niklasson -:cite:`t-Niklasson2008,t-Niklasson2014` is to introduce auxilliary electronic +:cite:`t-Niklasson2008,t-Niklasson2014` is to introduce auxiliary electronic degrees of freedom into the Lagrangian, which can be propagated via time-reversible integrators. @@ -138,7 +138,7 @@ model such systems by positing a set of equations of *non-Hamiltonian* equations of motion, and proving that they generate the correct statistical ensemble :cite:`t-Tuckerman2010`. This is the extended system approach: we modify the Hamiltonian to include the thermostat and/or barostat degrees of freedom, derive -the (pseudo-) Hamiltonian equations of motion, and demostrate that the correct +the (pseudo-) Hamiltonian equations of motion, and demonstrate that the correct phase space distribution for the ensemble is recovered. Go to :ref:`top `. @@ -192,7 +192,7 @@ expansion, This is recognisable as the velocity Verlet algorithm with extended Lagrangian integration which can be reduced to a single step, as described in -:ref:`Extended Lagrangian Born-Oppenheimer MD (XL-BOMD)`, with a half time step +:ref:`Extended Lagrangian Born-Oppenheimer MD (XL-BOMD) `, with a half time step integration of the Nose-Hoover chain equations of motion before and after. For full details of the integration scheme, see Hirakawa *et al*. :cite:`t-Hirakawa2017`. @@ -204,11 +204,11 @@ Go to :ref:`top `. Isobaric-Isothermal (NPT) ensemble ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -The Parinello-Rahman equations of motion :cite:`t-Parrinello1981` extend the +The Parrinello-Rahman equations of motion :cite:`t-Parrinello1981` extend the fixed cell equations of motion to include the cell degrees of freedom in the extended system approach. We use the Martyna-Tobias-Tuckerman-Klein modification :cite:`t-Martyna1996`, which couples the variable cell equations of motion to a -Nose-Hoover chain the themrostat the system, recovering the isobaric-isothermal +Nose-Hoover chain to thermostat the system, recovering the isobaric-isothermal (NPT) ensemble. For an unconstrained cell (i.e. the lattice vectors can change freely), the equations of motion are, @@ -218,7 +218,7 @@ freely), the equations of motion are, \mathbf{\dot{h}} &= \frac{\mathbf{p}_g\mathbf{h}}{W_g} \\ \mathbf{\dot{p}_g} &= V(\mathbf{P}_\mathrm{int}-\mathbf{I}P_\mathrm{ext}) + \left[\frac{1}{N_f}\sum_{i=1}^N\frac{\mathbf{p}_i^2}{m_i}\right]\mathbf{I} - \frac{p_\xi}{Q}\mathbf{p}_g \\ \dot{\xi} &= \frac{p_\xi}{Q} \\ - \mathbf{\dot{p}}_g &= \sum_{i=1}^N\frac{\mathbf{p}_i^2}{m_i} + \frac{1}{W_g}\mathrm{Tr}[\mathbf{p}_g^T\mathbf{p}_g] - (N_f + d^2)kT + \dot{p}_\xi &= \sum_{i=1}^N\frac{\mathbf{p}_i^2}{m_i} + \frac{1}{W_g}\mathrm{Tr}[\mathbf{p}_g^T\mathbf{p}_g] - (N_f + d^2)kT Here, :math:`\mathbf{r}_i`, :math:`\mathbf{p}_i` and :math:`m_i` are the position, momentum and mass of particle :math:`i` respectively, :math:`\xi`, @@ -281,17 +281,17 @@ Weak coupling thermostat/barostat ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Instead of modifying the Hamiltonian, the Berendsen-type weak coupling method -:cite:`t-Berendsen1984` involves coupling the ionic degrees of freedom to a an +:cite:`t-Berendsen1984` involves coupling the ionic degrees of freedom to an external temperature and/or pressure bath via "the principle of least local perturbation consistent with the required global coupling." Thermostatting is -acheived via a Langevin-type equation of motion, in which the system is globally -coupled to a heat bath and subjected to random noise: +achieved through deterministic global feedback from the instantaneous +temperature; unlike Langevin dynamics, no random force is applied: .. math:: m_i\ddot{\mathbf{r}}_i = \mathbf{F}_i + m_i \gamma\left(\frac{T_0}{T}-1\right)\dot{\mathbf{r}}_i, where :math:`\gamma` is a global friction constant chosen to be the same for all -particles. This can be acheived in practice by rescaling the velocities +particles. This can be achieved in practice by rescaling the velocities :math:`\mathbf{v}_i \rightarrow \lambda\mathbf{v}_i`, where :math:`\lambda` is, .. math:: @@ -306,13 +306,14 @@ be expressed, i.e. the fractional coordinates are scaled by a factor determined by the difference between the internal and external pressures, the isothermal -compressibility :math:`\beta` and a pressure coupling time constant $\tau_P$. +compressibility :math:`\beta` and a pressure coupling time constant +:math:`\tau_P`. In the isotropic case, the cell scaling factor :math:`\mu` can be expressed, .. math:: \mu = \left[ 1 - \frac{\Delta t}{\tau_P}(P_0 - P)\right]^{\frac{1}{3}}, -where the compressibility is absorbed into the time time constant +where the compressibility is absorbed into the time constant :math:`\tau_P`. Allowing for fluctuations of all cell degrees of freedom, the scaling factor becomes, @@ -350,7 +351,7 @@ process. Practically, the particle velocities are rescaled by a factor of Where :math:`R_i` is a set of :math:`N_f` normally distributed random numbers with unitary variance. This method can be applied to thermostat the NPT ensemble -by barostatting the system with the Parinello-Rahman method, and using the above +by barostatting the system with the Parrinello-Rahman method, and using the above expressions, but with additional :math:`R_i`'s for the cell degrees of freedom, and thermostatting the cell velocities as well as the particle velocities :cite:`t-Bussi2009`. diff --git a/src/DiagModule.f90 b/src/DiagModule.f90 index ca73905e4..6c2edba34 100644 --- a/src/DiagModule.f90 +++ b/src/DiagModule.f90 @@ -3171,12 +3171,13 @@ end subroutine occupy !! CREATION DATE !! 31/07/2026 !! MODIFICATION HISTORY - !! + !! 2026/09/15 16:03 dave + !! Added Fermi level output for linear scaling !! SOURCE subroutine write_gaps(spin_ch) use units - use global_module, only: nspin, flag_fix_spin_population + use global_module, only: nspin, flag_fix_spin_population, flag_diagonalisation, mu_DMM implicit none @@ -3196,7 +3197,7 @@ subroutine write_gaps(spin_ch) spin_st = 1 spin_end = 1 end if - if((flag_smear_type==0.or.flag_integer_occ).and.myid==0) then + if(flag_diagonalisation.and.(flag_smear_type==0.or.flag_integer_occ).and.myid==0) then do spin=spin_st,spin_end if(flag_gap(spin)) then if(nspin>1.and.flag_fix_spin_population) then!spin_end-spin_st>0) then @@ -3222,6 +3223,10 @@ subroutine write_gaps(spin_ch) end if write(io_lun,fmt='(4x,"Fermi level= ",f12.5," ",a2)') en_conv*Efermi(spin),en_units(energy_units) end do + else if(.not.flag_diagonalisation.and.myid==0) then + do spin=spin_st,spin_end + write(io_lun,fmt='(4x,"Fermi level= ",f12.5," ",a2)') en_conv*mu_DMM(spin),en_units(energy_units) + end do end if return end subroutine write_gaps diff --git a/src/UpdateInfo_module.f90 b/src/UpdateInfo_module.f90 index f835b50c2..d9cb50f14 100644 --- a/src/UpdateInfo_module.f90 +++ b/src/UpdateInfo_module.f90 @@ -1302,6 +1302,8 @@ end subroutine CommMat_send_data !! Added InfoGlob in dummy arguments, and removed glob2node_old !! 2026/02/10 10:45 dave !! Bug fix: changed nreq to single integer (used in block receive) + !! 2026/09/10 TM + !! Bug fix: changed the type of nreq !! SOURCE !! subroutine alloc_recv_array(InfoGlob,irecv_array,irecv2_array,recv_array, & @@ -1336,7 +1338,8 @@ subroutine alloc_recv_array(InfoGlob,irecv_array,irecv2_array,recv_array, & logical :: flag_find_old ! -- process 2. -- ! integer :: ind_nnd,ibeg,isize,inode_send - integer :: tag,ierr,nreq + integer :: tag,ierr + integer :: nreq(MPI_STATUS_SIZE) integer, allocatable :: isort_node(:),irecv_array(:) ! -- process 3. -- ! integer :: iprim_remote,isize1,isize2,ia,nalpha,nj,njbeta,iglob_local diff --git a/src/initial_read_module.f90 b/src/initial_read_module.f90 index 83ef0b3d6..6e08cd63b 100644 --- a/src/initial_read_module.f90 +++ b/src/initial_read_module.f90 @@ -13,7 +13,7 @@ !! Holds the initial read routines !! USES !! atoms, construct_module, datatypes, DiagModule, dimens, fdf, GenComms, global_module, -!! group_module, io_module, maxima_module, numbers, parse, primary_module, +!! group_module, io_module, maxima_module, numbers, parse, primary_module, !! pseudopotential_data, species_module !! AUTHOR !! D.R.Bowler @@ -24,11 +24,11 @@ !! Bug fix (from TM) for read_and_write !! 17/06/2002 dave !! Added read to find solution method - order N or exact diagonalisation -!! 15:55, 25/09/2002 mjg & drb +!! 15:55, 25/09/2002 mjg & drb !! Added flag for initial charge from initial K -!! 10:41, 06/03/2003 drb +!! 10:41, 06/03/2003 drb !! Added tags for linear mixing for self-consistency in read_input -!! 11:02, 24/03/2003 drb +!! 11:02, 24/03/2003 drb !! Simplified read_and_write !! 14:51, 2003/06/09 dave !! Added new way of reading atomic coordinates, independent of make_prt @@ -53,9 +53,9 @@ !! 2015/06/11 16:30 lat !! Added possibility of user to define filename for PAO/pseudo *ion !! Fixed EXX variable definitions -!! Try to make blocks of variables clearer +!! Try to make blocks of variables clearer !! Added experimental backtrace -!! +!! !! SOURCE !! module initial_read @@ -76,13 +76,13 @@ module initial_read ! ------------------------------------------------------------------------------ ! Subroutine read_and_write ! ------------------------------------------------------------------------------ - + !!****f* initial_read/read_and_write * !! - !! NAME + !! NAME !! read_and_write !! USAGE - !! + !! !! PURPOSE !! Reads the input data and writes out the parameters for the run !! INPUTS @@ -108,11 +108,11 @@ module initial_read !! 31/05/2002 dave !! Bug fix (TM) - changed np to ind_part in loop to build atom !! positions from make_prt.dat - !! 11:02, 24/03/2003 drb + !! 11:02, 24/03/2003 drb !! Removed variables associated with old phi initialisation !! 14:52, 2003/06/09 dave !! Removed atomic coordinate fiddling - !! 14:59, 02/05/2005 dave + !! 14:59, 02/05/2005 dave !! Added calculation of number of electrons in cell !! 2007/05/08 17:00 dave !! Added net charge on cell to number_of_bands for consistency @@ -128,7 +128,7 @@ module initial_read !! - Changed spin implementation !! 2014/02/04 M.Arita !! - Added call for initial_read_aux for constraint-MD - !! 2015/05/11 lat + !! 2015/05/11 lat !! - Added optional total spin magnetization !! 2015/06/08 lat !! - Added experimental backtrace @@ -195,8 +195,8 @@ subroutine read_and_write(start, start_L, inode, ionode, & use read_gto_info, only: read_gto_new !use gto_format, only: gto use gto_format_new, only: gto - - use cdft_data, only: cDFT_NAtoms, & + + use cdft_data, only: cDFT_NAtoms, & cDFT_NumberAtomGroups, cDFT_AtomList use memory_module, only: reg_alloc_mem, type_dbl use primary_module, only: bundle, make_prim @@ -291,7 +291,7 @@ subroutine read_and_write(start, start_L, inode, ionode, & if(flag_functional_type/=pseudo(i_species)%functional) then if(flag_different_functional) then call cq_warn(sub_name, "Functional in input file differs to pseudopotential but proceeding: ",& - flag_functional_type, pseudo(i_species)%functional) + flag_functional_type, pseudo(i_species)%functional) else call cq_abort("Functional in input file differs to pseudopotential: ",& flag_functional_type, pseudo(i_species)%functional) @@ -308,7 +308,7 @@ subroutine read_and_write(start, start_L, inode, ionode, & ! allocate(gto_file(n_species),STAT=stat) if(stat /= 0) call cq_abort("Error allocating gto_file in read_and_write: ",n_species,stat) - allocate(gto(n_species),STAT=stat) + allocate(gto(n_species),STAT=stat) if(stat /= 0) call cq_abort ("Error allocating gto in read_and_write",stat) ! call read_gto_new(inode,ionode,n_species) @@ -442,8 +442,8 @@ subroutine read_and_write(start, start_L, inode, ionode, & ! multi-site SFs (MSSFs) are symmetry-breaking usually. ! The default of the number of the MSSFs is SingleZeta-size. ! SZP-size MSSFs with flag_MSSF_nonminimal is also available. - ! If the user wants to use the other number of MSSFs, - ! the user must provide the initial SFcoeffmatrix or the trial vectors for the LFD method. + ! If the user wants to use the other number of MSSFs, + ! the user must provide the initial SFcoeffmatrix or the trial vectors for the LFD method. if (nsf_species(i).eq.count_SZ) then MSSF_nonminimal_species(i) = 1 ! SZ-size MSSF else @@ -473,7 +473,7 @@ subroutine read_and_write(start, start_L, inode, ionode, & endif endif endif - else + else ! If number of support functions is less than total number of ang. mom. components (ignoring ! for now multiple zetas) then there is a formal problem with basis set: we require the user ! to set an additional flag to assert that this is really desired @@ -529,7 +529,7 @@ subroutine read_and_write(start, start_L, inode, ionode, & if(flag_plusUproj_atom(i_species)) then do i_ang = 0, pao(i_species)%greatest_angmom if(pao(i_species)%angmom(i_ang)%n_zeta_in_angmom>0) then - if (i_ang.eq.l_plusUproj) then + if (i_ang.eq.l_plusUproj) then ! set U parameters for PAO with (n_plusUproj,l_plusUproj,z_plusUproj) z_temp = 0 do i_zeta = 1, pao(i_species)%angmom(i_ang)%n_zeta_in_angmom @@ -606,7 +606,7 @@ subroutine read_and_write(start, start_L, inode, ionode, & if(flag_InitialAtomicSpin) then do i = 1, n_species charge_tmp = charge_up(i) + charge_dn(i) - if(charge_tmp < RD_ERR) then ! + if(charge_tmp < RD_ERR) then ! charge_up(i) = half*charge(i) charge_dn(i) = half*charge(i) !if(inode .eq. ionode) write(io_lun,fmt='(6x,a,i3,a,2f15.8)') & @@ -621,12 +621,12 @@ subroutine read_and_write(start, start_L, inode, ionode, & end do !if(inode .eq. ionode) write(io_lun,fmt='(6x,a,2f15.8)') 'ne_spin_in_cell(1:2) = ',ne_spin_in_cell(1),ne_spin_in_cell(2) endif - ! + ! ! ! Set up various lengths, volumes, reciprocals etc. for convenient use call set_dimensions(inode, ionode,HNL_fac, non_local, n_species, & non_local_species, core_radius) - ! + ! ! write out some information on the run if (inode == ionode) & call write_info(titles, mu, vary_mu, HNL_fac, numprocs) @@ -648,18 +648,18 @@ end subroutine read_and_write ! ------------------------------------------------------------------------------ ! Subroutine read_input ! ------------------------------------------------------------------------------ - + !!****f* initial_read/read_input * !! - !! NAME + !! NAME !! read_input !! USAGE - !! + !! !! PURPOSE !! reads the input file !! INPUTS - !! - !! + !! + !! !! USES !! datatypes, global_module, atoms, dimens, species_module, !! pseudopotential_data, GenComms, fdf, parse @@ -680,27 +680,27 @@ end subroutine read_and_write !! Incorporated into initial_read module !! 29/05/2002 dave !! Added scan for SolutionMethod tag - !! 15:54, 25/09/2002 mjg & drb + !! 15:54, 25/09/2002 mjg & drb !! Added scan for make_initial_charge_from_K (which sets find_chdens) !! 16:31, 2003/02/03 dave !! Added scans for type of run (e.g. MD, CG etc) !! 07:49, 2003/02/04 dave !! Added boolean variables to control run (vary blips, self-consistency) - !! 10:41, 06/03/2003 drb + !! 10:41, 06/03/2003 drb !! Added tags for linear mixing in self-consistency - !! 10:09, 12/03/2003 drb + !! 10:09, 12/03/2003 drb !! Added an end tag !! 14:52, 2003/06/09 dave !! Added new way to read atomic coordinates and preconditioning flag !! 08:31, 2003/10/01 dave !! Changed flag_vary_blips to flag_vary_basis - !! 13:15, 03/10/2003 drb + !! 13:15, 03/10/2003 drb !! Bug fix - missing . !! 2004/10/29 drb !! Added various new flags for self consistency - !! 14:59, 02/05/2005 dave + !! 14:59, 02/05/2005 dave !! Added call to check for net charge - initialises ne_in_cell - !! 09:14, 11/05/2005 dave + !! 09:14, 11/05/2005 dave !! Store max_L_iterations in global variable !! 13:34, 2006/07/09 ast !! Added flag for functional @@ -712,22 +712,22 @@ end subroutine read_and_write !! Added keyword MaxEfIter to stop infinite loop in findFermi !! 2008/01/24 Veronika !! Changed the default of General.ManyProcessors to .true. - !! 12:18, 14/02/2008 drb + !! 12:18, 14/02/2008 drb !! Added options for buffer around primary and covering sets !! 2008/07/16 ast !! New keywords for timers !! 2009/07/08 16:47 dave !! Added keyword for one-to-one PAO to SF assigment !! 2009/07/23 ast - !! More general routine to name timer files, so that they are not + !! More general routine to name timer files, so that they are not !! to a maximum number of processes, as it was the case before, !! with only three figures (max. 999 processes) !! 2010/06/18 lt !! Added control flags for Methfessel-Paxton smearing method !! 2010/07/22 15.53 Lianheng - !! Moved control flags associated with Diagonalisation method from - !! here to readDiagInfo() subroutine. Added Methfessel-Paxton - !! smearing and related control flags. Changed SC.MaxEfIter to + !! Moved control flags associated with Diagonalisation method from + !! here to readDiagInfo() subroutine. Added Methfessel-Paxton + !! smearing and related control flags. Changed SC.MaxEfIter to !! Diag.MaxEfIter (moved to readDiagInfo()) for greater consistency !! 2011/07/21 16:50 dave !! Added input flags for cDFT @@ -861,7 +861,7 @@ end subroutine read_and_write !! AtomMove.ReuseL => AtomMove.ReuseLorK => AtomMove.ReuseDM !! 2020/01/06 15:43 dave !! Keywords for equilibration - !! 2020/01/07 tsuyoshi + !! 2020/01/07 tsuyoshi !! Default setting of MakeInitialChargeFromK has been changed !! 2020/12/14 lionel !! EXX: added filtering option for EXX and cleaning @@ -874,7 +874,7 @@ end subroutine read_and_write !! 2022/12/14 10:01 dave and tsuyoshi !! Update test for solution method (diagon vs ordern) following issue #47 !! 2024/05/22 14:40 nakata - !! Added DFT+U related input and output + !! Added DFT+U related input and output !! 2024/12/03 lionel !! Added grid specification of EXX coarse/standard/fine !! 2025/02/03 nakata @@ -940,7 +940,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& flag_write_projected_DOS, & E_wf_min, E_wf_max, flag_wf_range_Ef, & mx_temp_matrices, flag_neutral_atom, flag_diagonalisation, & - flag_DFTplusU, flag_first_diag, flag_write_occ_mat, & + flag_DFTplusU, flag_first_diag, flag_write_occ_mat, & flag_SpinDependentSF, flag_Multisite, flag_LFD, flag_SFcoeffReuse, & flag_opt_cell, cell_constraint_flag, flag_variable_cell, & cell_en_tol, optcell_method, cell_stress_tol, & @@ -1128,19 +1128,19 @@ subroutine read_input(start, start_L, titles, vary_mu,& else io_lun = 6 end if - ! + ! if (fdf_boolean('IO.WriteOutToASEFile',.false.)) then !call io_assign(io_ase) ! Reserve this unit on all processes !write_ase = .false. if (inode == ionode) then write_ase = .true. - call io_assign(io_ase) + call io_assign(io_ase) open(unit=io_ase,file=ase_file,iostat=stat) if (stat /= 0) & call cq_abort("Failed to open ASE Conquest output file", stat) end if else - if (inode == ionode) write_ase = .false. + if (inode == ionode) write_ase = .false. end if !!$ !!$ @@ -1149,7 +1149,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& !!$ !!$ !!$ - BackTraceOn = fdf_boolean('IO.BackTraceOn',.false.) + BackTraceOn = fdf_boolean('IO.BackTraceOn',.false.) ! Where will the timers will be written? TimingOn = fdf_boolean('IO.TimingOn',.false.) ! tmr_rmv001 if(TimingOn) then ! tmr_rmv001 @@ -1176,10 +1176,10 @@ subroutine read_input(start, start_L, titles, vary_mu,& end if ! ionode switch ! tmr_rmv001 end if ! tmr_rmv001 end if ! tmr_rmv001 -!!$ -!!$ **** +!!$ +!!$ **** !!$ I think we may need (std) output file to be open before this point -!!$ at the very beginning in the main.f90 +!!$ at the very beginning in the main.f90 !!$ !!$ !!$ @@ -1252,7 +1252,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& restart_T = fdf_boolean('General.LoadInvS',.false.) ! Is there a net charge on the cell ? ne_in_cell = fdf_double('General.NetCharge',zero) - ! Read coordinates file + ! Read coordinates file flag_fractional_atomic_coords = & fdf_boolean('IO.FractionalAtomicCoords',.true.) flag_check_init_atomic_coords = & @@ -1267,7 +1267,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& ! ! Spin polarised calculation ! - if we are doing spin polarised calculations or not - flag_spin_polarisation = fdf_boolean('Spin.SpinPolarised', .false.) + flag_spin_polarisation = fdf_boolean('Spin.SpinPolarised', .false.) flag_fix_spin_population = fdf_boolean('Spin.FixSpin', .false.) if (flag_spin_polarisation) then nspin = 2 @@ -1321,13 +1321,13 @@ subroutine read_input(start, start_L, titles, vary_mu,& in_block_y = fdf_integer('Grid.InBlockY',4) in_block_z = fdf_integer('Grid.InBlockZ',4) ! Solution method - O(N) or diagonalisation - diagonalisation default - method = fdf_string(6,'DM.SolutionMethod','diagon') + method = fdf_string(6,'DM.SolutionMethod','diagon') if(leqi(method,'ordern')) then flag_diagonalisation = .false. flag_check_Diag = .false. else flag_diagonalisation = .true. - flag_check_Diag = .true. + flag_check_Diag = .true. end if ! Read basis set basis_string = fdf_string(10,'Basis.BasisSet','PAOs') @@ -1357,7 +1357,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& !!$ !!$ ! Hamann is default to fit with Conquest ion file generator - ps_type = fdf_string(5,'General.PseudopotentialType','haman') + ps_type = fdf_string(5,'General.PseudopotentialType','haman') ! Write out pseudopotential type if(leqi(ps_type,'siest')) then !if(inode==ionode.AND.iprint_init>0) & @@ -1377,7 +1377,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& inode==ionode) then write(io_lun,& fmt='(10x,"Setting FlagNewAngular to T for Siesta/& - &Abinit pseudopotentials")') + &Abinit pseudopotentials")') flag_angular_new = .true. end if ! Should we use neutral atom potential ? @@ -1412,7 +1412,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& call cq_warn(sub_name,"DFT+U not available with blips: turning off DFT+U") flag_DFTplusU = .false. end if - + !!! nakata DFT+U end !!$ !!$ @@ -1434,7 +1434,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& species_from_files = fdf_boolean('General.PAOFromFiles',.false.) flag_ghost=.false. if(fdf_block('ChemicalSpeciesLabel')) then - if(1+block_end-block_start RD_ERR) then @@ -1589,7 +1589,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& ! At present, we cannot use neutral atom potential for ghost atoms. ! 2018/Sep/26 Tsuyoshi Miyazaki ! Fixed for Neutral Atom Potential 2018/Nov/16 TM - ! but not for neutral atom projectors ... + ! but not for neutral atom projectors ... if(flag_neutral_atom_projector) then if(inode==ionode) write(io_lun,'(10x,A)') & 'Now, we cannot use neutral atom projector for ghost atoms. Please set General.NeutralAtomProjector as false.' @@ -1612,11 +1612,11 @@ subroutine read_input(start, start_L, titles, vary_mu,& ! enddo !endif endif -!!$ -!!$ -!!$ -!!$ -!!$ +!!$ +!!$ +!!$ +!!$ +!!$ ! Set variables for multisite support functions if (flag_Multisite) then flag_MSSF_nonminimal = fdf_boolean('Multisite.nonminimal', .false.) !nonmin_mssf @@ -1687,7 +1687,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& !!$ !!$ !!$ -!!$ +!!$ !blip_width = support_grid_spacing * ! fdf_double('blip_width_over_support_grid_spacing',four) if (leqi(runtype,'md')) then @@ -1744,7 +1744,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& flag_vary_basis = fdf_boolean('minE.VaryBasis', .false.) if(.NOT.flag_vary_basis) then flag_precondition_blips = .false. - else + else flag_precondition_blips = fdf_boolean('minE.PreconditionBlips',.true.) end if InitStep_paomin = fdf_double ('minE.InitStep_paomin', 5.0_double) @@ -1753,7 +1753,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& ! DRB 2018/02/26 turn off mixed L-SCF with diagonalisation if(flag_mix_L_SC_min .and. flag_diagonalisation) then flag_mix_L_SC_min = .false. - else if(flag_mix_L_SC_min .and. flag_self_consistent) then + else if(flag_mix_L_SC_min .and. flag_self_consistent) then flag_self_consistent = .false. end if ! Tweak 2007/03/23 DRB Make Pulay mixing default @@ -1776,8 +1776,8 @@ subroutine read_input(start, start_L, titles, vary_mu,& ! New residual flags jtlp 08/2019 flag_newresidual = fdf_boolean('SC.AbsResidual', .false.) flag_newresid_abs = fdf_boolean('SC.AbsResidual.Fractional', .true.) - ! When constructing charge density from K matrix at last step, we check the total - ! number of electrons. If the error of electron number (per total electron number) + ! When constructing charge density from K matrix at last step, we check the total + ! number of electrons. If the error of electron number (per total electron number) ! is larger than the following value, we use atomic charge density. (in update_H) threshold_resetCD = fdf_double('SC.Threshold.Reset',0.1_double) ! Surface dipole correction parameters @@ -1931,7 +1931,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& end if !!$ !!$ -!!$ D e l t a - S C F +!!$ D e l t a - S C F !!$ !!$ ! DeltaSCF flags @@ -2128,13 +2128,13 @@ subroutine read_input(start, start_L, titles, vary_mu,& !!$ !!$ !!$ -!!$ +!!$ flag_functional_type = fdf_integer('General.FunctionalType', 0) ! Read from pseudopotentials flag_different_functional = fdf_boolean('General.DifferentFunctional',.false.) ! Use different functional to ion files !!$ !!$ -!!$ E X X -!!$ +!!$ E X X +!!$ !!$ if ( flag_functional_type == functional_hyb_pbe0 ) then flag_exx = .true. @@ -2149,12 +2149,12 @@ subroutine read_input(start, start_L, titles, vary_mu,& r_exx = fdf_double ('EXX.Xrange', zero) ! else - ! don't touch we need it because matX is setup in set_dimensions + ! don't touch we need it because matX is setup in set_dimensions ! whatever is flag_exx flag_exx = .false. exx_scf = fdf_integer('EXX.MethodSCF', -1) - r_exx = one - ! + r_exx = one + ! end if ! if ( flag_exx ) then @@ -2179,7 +2179,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& exx_grid = fdf_string (20,'EXX.Grid','standard') exx_hgrid = fdf_double ('EXX.GridSpacing',zero) exx_radius = fdf_double ('EXX.IntegRadius',zero) - exx_scheme = fdf_integer('EXX.Scheme', 1 ) + exx_scheme = fdf_integer('EXX.Scheme', 1 ) exx_debug = fdf_boolean('EXX.Debug', .false. ) exx_overlap= fdf_boolean('EXX.Overlap',.true. ) ! @@ -2196,7 +2196,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& exx_pscheme = fdf_string (20,'EXX.FFTWSolver','ewald') if(exx_pscheme == 'ewald') then ewald_alpha = fdf_double('EXX.FFTWEwaldAlpha',3.0_double) - end if + end if end if !exx_mem = 1 end if @@ -2277,7 +2277,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& ! !!$ !!$ -!!$ +!!$ !!$ !!$ append_coords = fdf_boolean('AtomMove.AppendCoords',.true.) @@ -2327,7 +2327,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& ! DensityMatrix = 0; AtomicCharge = 1; LastStep = 2 method_UpdateChargeDensity = fdf_integer('AtomMove.InitialChargeDensity',DensityMatrix) - ! The keywords ( SC.ResetDensOnAtomMove and Multisite.LFD.UpdateWithAtomicDensity ) + ! The keywords ( SC.ResetDensOnAtomMove and Multisite.LFD.UpdateWithAtomicDensity ) ! should be removed in the near future. flag_reset_dens_on_atom_move = fdf_boolean('SC.ResetDensOnAtomMove',.false.) if(flag_reset_dens_on_atom_move) then @@ -2357,7 +2357,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& flag_XLBOMD = fdf_boolean('AtomMove.ExtendedLagrangian',.false.) ! zamaan 2018/03/03 I can't imagine a case where you would want to - ! restart a MD run without loading the various matrices, so I'm + ! restart a MD run without loading the various matrices, so I'm ! Defaulting some flags to true. Calling fdf to ensure that input.log ! remains consistent if (flag_MDcontinue) then @@ -2369,7 +2369,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& else flag_read_velocity = fdf_boolean('AtomMove.ReadVelocity',.false.) restart_DM = fdf_boolean('General.LoadDM', .false.) - if(restart_DM) then + if(restart_DM) then find_chdens = fdf_boolean('SC.MakeInitialChargeFromK',.true.) ! Allow DFT+U Hamiltonian with DM we have just loaded if(flag_DFTplusU) flag_first_diag = .false. @@ -2448,7 +2448,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& &positive integer:', constraints%n_grp) allocate (constraints%grp_name(constraints%n_grp) , & constraints%n_atom_in_grp(constraints%n_grp) , & - constraints%n_subgrp(constraints%n_grp)) + constraints%n_subgrp(constraints%n_grp)) do i = 1, constraints%n_grp if (fdf_block('RigidBondsGroups')) then read (unit=input_array(block_start+i-1),fmt=*) & @@ -2471,7 +2471,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& case('nve') md_thermo_type = fdf_string(20, 'MD.Thermostat', 'none') md_baro_type = fdf_string(20, 'MD.Barostat', 'none') - case('nvt') + case('nvt') md_thermo_type = fdf_string(20, 'MD.Thermostat', 'nhc') md_baro_type = fdf_string(20, 'MD.Barostat', 'none') case('npt') @@ -2481,14 +2481,14 @@ subroutine read_input(start, start_L, titles, vary_mu,& end select md_tau_T = fdf_double('MD.tauT', -one) md_tau_T_equil = fdf_double('MD.tauTEquil', one) - md_n_nhc = fdf_integer('MD.nNHC', 5) + md_n_nhc = fdf_integer('MD.nNHC', 5) md_n_ys = fdf_integer('MD.nYoshida', 1) md_n_mts = fdf_integer('MD.nMTS', 1) flag_thermoDebug = fdf_boolean('MD.ThermoDebug',.false.) md_t_drag = fdf_double('MD.TDrag', zero) if (leqi(md_thermo_type, 'nhc')) then - allocate(md_nhc_mass(md_n_nhc)) - allocate(md_nhc_cell_mass(md_n_nhc)) + allocate(md_nhc_mass(md_n_nhc)) + allocate(md_nhc_cell_mass(md_n_nhc)) md_nhc_mass = one md_nhc_cell_mass = one if (fdf_block('MD.NHCMass')) then @@ -2595,7 +2595,7 @@ end subroutine read_input !! AUTHOR !! T. Miyazaki !! CREATION DATE - !! 2017/11/03 + !! 2017/11/03 !! MODIFICATION HISTORY !! 2020/07/30 tsuyoshi !! - Moved from read_input @@ -2605,10 +2605,10 @@ end subroutine read_input ! this subroutine checks the compatibility between keywords defined in read_input ! 2017.11(Nov).03 Tsuyoshi Miyazaki - ! + ! ! we don't need to worry about which parameter is defined first. ! - subroutine check_compatibility + subroutine check_compatibility use global_module, only: flag_move_atom, ni_in_cell, & runtype, flag_XLBOMD, flag_diagonalisation, & flag_LmatrixReuse, flag_SFcoeffReuse, flag_DumpMatrices, & @@ -2616,7 +2616,7 @@ subroutine check_compatibility use input_module, only: leqi use density_module,only: method_UpdateChargeDensity,DensityMatrix,AtomicCharge use GenComms, only: cq_warn - + implicit none logical :: flag_FixedAtoms integer :: ig, k @@ -2643,7 +2643,7 @@ subroutine check_compatibility 'AtomMove.InitialChargeDensity is changed to AtomicCharge, since AtomMove.ReuseDM is false') endif - !flag_DumpMatrices : at present, we need matrix files to reuse the previous matrix data + !flag_DumpMatrices : at present, we need matrix files to reuse the previous matrix data if(.not.flag_DumpMatrices) then if(flag_SFcoeffReuse .or. flag_LmatrixReuse) then flag_DumpMatrices = .true. @@ -2674,20 +2674,20 @@ end subroutine check_compatibility ! ------------------------------------------------------------------------------ ! Subroutine allocate_species_vars ! ------------------------------------------------------------------------------ - + !!****f* initial_read/allocate_species_vars * !! - !! NAME + !! NAME !! allocate_species_vars !! USAGE - !! + !! !! PURPOSE !! Allocates variables which depend on number of species !! INPUTS - !! - !! + !! + !! !! USES - !! + !! !! AUTHOR !! D. R. Bowler !! CREATION DATE @@ -2718,7 +2718,7 @@ subroutine allocate_species_vars use memory_module, only: reg_alloc_mem, type_dbl use species_module, only: n_species, nsf_species, nlpf_species, npao_species, natomf_species, & charge, charge_up, charge_dn - use species_module, only: mass, non_local_species, ps_file, ch_file, phi_file + use species_module, only: mass, non_local_species, ps_file, ch_file, phi_file use species_module, only: species_label, species_file, type_species use global_module, only: area_general use GenComms, only: cq_abort @@ -2817,13 +2817,13 @@ end subroutine allocate_species_vars ! ------------------------------------------------------------------------------ ! Subroutine write_info ! ------------------------------------------------------------------------------ - + !!****f* initial_read/write_info * !! - !! NAME + !! NAME !! write_info !! USAGE - !! + !! !! PURPOSE !! Writes out information about the run !! INPUTS @@ -2835,7 +2835,7 @@ end subroutine allocate_species_vars !! CREATION DATE !! 9/11/93 (!) !! MODIFICATION HISTORY - !! 10/8/01 DRB + !! 10/8/01 DRB !! Changed to F90 format and tidied !! 10/05/2002 dave !! Changed titles to be a single line @@ -2966,10 +2966,10 @@ subroutine write_info(titles, mu, vary_mu, HNL_fac, NODES) end if ! Ground state search details write(io_lun,fmt='(/4x,"Ground state search:")') - if(flag_basis_set==blips) then + if(flag_basis_set==blips) then write(io_lun,'(6x,"Support functions represented with blip basis")') else - write(io_lun,'(6x,"Support functions represented with PAO basis")') + write(io_lun,'(6x,"Support functions represented with PAO basis")') if(flag_Multisite) then if(flag_LFD) then write(io_lun,'(6x,"Multi-site SFs used with local filter diagonalisation")') @@ -3006,7 +3006,7 @@ subroutine write_info(titles, mu, vary_mu, HNL_fac, NODES) end select end if else - write(io_lun,fmt='(6x,"Solving for the K matrix using ",a16)') 'order N with LNV' + write(io_lun,fmt='(6x,"Solving for the K matrix using ",a16)') 'order N with LNV' end if if(iprint_init>0) write(io_lun,fmt='(/4x,"Integration grid size: ",i4," x ",i4," x ",i4)') & n_grid_x, n_grid_y, n_grid_z @@ -3015,10 +3015,10 @@ subroutine write_info(titles, mu, vary_mu, HNL_fac, NODES) in_block_x, in_block_y, in_block_z write(io_lun,fmt='(/4x,"Integration grid spacing: ",f6.3,a3," x",f6.3,a3," x",f6.3,a3)') & - dist_conv*(r_super_x/n_grid_x), d_units(dist_units), & + dist_conv*(r_super_x/n_grid_x), d_units(dist_units), & dist_conv*(r_super_y/n_grid_y), d_units(dist_units), & dist_conv*(r_super_z/n_grid_z), d_units(dist_units) - + ! print in Conquest output write(io_lun,fmt='(/4x,"Number of species: ",i2)') n_species write(io_lun,fmt='(4x,a56)') '--------------------------------------------------------' @@ -3031,7 +3031,7 @@ subroutine write_info(titles, mu, vary_mu, HNL_fac, NODES) n, mass(n), charge(n), dist_conv*RadiusSupport(n), nsf_species(n), species_label(n) end do write(io_lun,fmt='(4x,a56)') '--------------------------------------------------------' - + ! ! BEGIN %%%% ASE printing %%%% ! @@ -3044,7 +3044,7 @@ subroutine write_info(titles, mu, vary_mu, HNL_fac, NODES) write(io_ase,fmt='(4x,"| # mass (au) Charge (e) SF Rad (",a2,") NSF Label |")') & d_units(dist_units) write(io_ase,fmt='(4x,a56)') '--------------------------------------------------------' - + do n=1, n_species write(io_ase,fmt='(4x,"|",i4,2x,f9.3,3x,f9.3,4x,f9.3,2x,i3,2x,a7,"|")') & n, mass(n), charge(n), dist_conv*RadiusSupport(n), nsf_species(n), species_label(n) @@ -3061,15 +3061,15 @@ subroutine write_info(titles, mu, vary_mu, HNL_fac, NODES) if(iprint_init>0) then if(flag_vary_basis) then - write(io_lun,fmt='(10x,"Support function tolerance: ",f12.8)') energy_tolerance + write(io_lun,fmt='(10x,"Support function tolerance: ",e10.3)') energy_tolerance write(io_lun,fmt='(10x,"Support function iterations: ",i4)') n_support_iterations end if if(flag_self_consistent) then - write(io_lun,fmt='(10x,"SCF tolerance: ",f12.8)') sc_tolerance + write(io_lun,fmt='(10x,"SCF tolerance: ",e10.3)') sc_tolerance write(io_lun,fmt='(10x,"SCF iterations: ",i4)') maxitersSC end if if(.not.flag_diagonalisation) then - write(io_lun,fmt='(10x,"O(N) tolerance: ",f12.8)') L_tolerance + write(io_lun,fmt='(10x,"O(N) tolerance: ",e10.3)') L_tolerance write(io_lun,fmt='(10x,"O(N) iterations: ",i4)') n_L_iterations end if end if @@ -3122,7 +3122,7 @@ subroutine write_info(titles, mu, vary_mu, HNL_fac, NODES) !!! 2024.05.20 nakata DFT+U if(flag_DFTplusU) then - write(io_lun,fmt='(10x,"DFT+U calculation will be performed")') + write(io_lun,fmt='(10x,"DFT+U calculation will be performed")') write(io_lun,fmt='(/4x,"Number of DFT-U projectors: ",i2)') num_plusUproj @@ -3149,10 +3149,10 @@ end subroutine write_info ! ----------------------------------------------------------------------------- ! Subroutine readDiagInfo ! ----------------------------------------------------------------------------- - + !!****f* DiagModule/readDiagInfo * !! - !! NAME + !! NAME !! readDiagInfo !! USAGE !! readDiagInfo(Number of Kpts, proc grid rows, proc grid cols, @@ -3191,7 +3191,7 @@ end subroutine write_info !! must contain the same number of processes. !! 2011/09/02 L.Tong !! Corrected Conquest default values for ScaLAPACK processor grid - !! definitions when k-point parallelisation is on. + !! definitions when k-point parallelisation is on. !! Added local variable proc_per_group !! 2011/09/05 L.Tong !! Added print out input information for k-point parallelisation @@ -3257,7 +3257,7 @@ subroutine readDiagInfo real(double) :: a, sum, dkx, dky, dkz, dk integer :: proc_per_group logical :: ms_is_prime - + ! k-point mesh type logical :: mp_mesh, done, flag_lines_kpoints, flag_gamma, test_ase integer, dimension(1:3) :: mp @@ -3268,7 +3268,7 @@ subroutine readDiagInfo integer :: counter character(len=2) :: suffix - !****lat<$ + !****lat<$ call start_backtrace(t=backtrace_timer,who='readDiagInfo',where=1,level=2) !****lat>$ @@ -3301,7 +3301,7 @@ subroutine readDiagInfo proc_groups = 1 end if - if (iprint_init > 1 .AND. inode==ionode) write (io_lun, 11) proc_groups + if (iprint_init > 1 .AND. inode==ionode) write (io_lun, 11) proc_groups ! Read/choose ScaLAPACK processor grid dimensions if(fdf_defined('Diag.ProcRows')) then proc_rows = fdf_integer('Diag.ProcRows',0) @@ -3314,7 +3314,7 @@ subroutine readDiagInfo proc_per_group = numprocs / proc_groups if(proc_rows*proc_cols > proc_per_group) & call cq_abort('Diag: proc grid product is too large: ', & - proc_rows*proc_cols, proc_per_group) + proc_rows*proc_cols, proc_per_group) else ! default values if(numprocs<9) then ! Recommended by ScaLapack User Guide @@ -3327,7 +3327,7 @@ subroutine readDiagInfo proc_per_group = numprocs / proc_groups a = sqrt(real(proc_per_group, double)) proc_rows = floor(a)+1 - do while(.NOT.done.AND.proc_rows>1) + do while(.NOT.done.AND.proc_rows>1) proc_rows = proc_rows - 1 proc_cols = floor(real(proc_per_group,double) / & real(proc_rows,double)) @@ -3341,12 +3341,12 @@ subroutine readDiagInfo do i=1,ni_in_cell matrix_size = matrix_size + nsf_species(species_glob(i)) end do - + ! Test if matrix_size is a prime number ms_is_prime = is_prime(matrix_size) if ( ms_is_prime ) call cq_warn(sub_name,'matrix size is a prime number', matrix_size) - - ! padH or not :temporary? + + ! padH or not :temporary? flag_padH = fdf_boolean('Diag.PaddingHmatrix',.true.) if(flag_padH) then @@ -3363,7 +3363,7 @@ subroutine readDiagInfo flag_padH = .false. endif endif - + if(.not.flag_padH) then if(fdf_defined('Diag.BlockSizeR')) then block_size_r = fdf_integer('Diag.BlockSizeR',1) @@ -3381,10 +3381,10 @@ subroutine readDiagInfo block_size_c = block_size_r call cq_warn(sub_name,'Use of PaddingHmatrix is recommended.') call cq_warn(sub_name,'prime: set block_size_c = block_size_r = 1 ') - else ! + else ! done = .false. block_size_r = matrix_size/max(proc_rows,proc_cols)+1 - do while(.NOT.done) + do while(.NOT.done) block_size_r = block_size_r - 1 a = real(matrix_size)/real(block_size_r) if(a - real(floor(a))<1e-8_double) done = .true. @@ -3404,7 +3404,7 @@ subroutine readDiagInfo write(io_lun,2) block_size_r, block_size_c write(io_lun,3) proc_rows, proc_cols end if - + !Using ELPA or not flag_use_elpa = fdf_boolean('Diag.UseELPA',.false.) @@ -3573,7 +3573,7 @@ subroutine readDiagInfo end if else ! Read Monkhorst-Pack mesh coefficients - ! Default is Gamma point only + ! Default is Gamma point only if(iprint_init>1.AND.inode==ionode) & write(io_lun,fmt='(/8x,"Using Monkhorst-Pack Kpoint mesh"//)') flag_gamma = fdf_boolean('Diag.GammaCentred',.false.) @@ -3645,7 +3645,7 @@ subroutine readDiagInfo wtk_tmp(1:nkp_tmp) = one ! Generate fractional k-point coordinates plus shift ! Assume orthorhombic cell for now - do i = 1, mp(1) ! x axis + do i = 1, mp(1) ! x axis do j = 1, mp(2) ! y axis do k = 1, mp(3) ! z axis counter = k + (j-1) * mp(3) + (i-1) * mp(2) * mp(3) @@ -3683,13 +3683,13 @@ subroutine readDiagInfo end if end do - ! How many k-points are now left? + ! How many k-points are now left? nkp = 0 - ! DRB fix + ! DRB fix sum = zero do i = 1, nkp_tmp if (wtk_tmp(i) > RD_ERR ) then - nkp = nkp + 1 + nkp = nkp + 1 sum = sum + wtk_tmp(i) end if end do @@ -3726,7 +3726,7 @@ subroutine readDiagInfo else if(iprint_init>1.AND.inode==ionode) then write(io_lun,fmt='(8x,i4,a)') nkp, ' symmetry inequivalent Kpoints' end if - + do i = 1, nkp kk(1,i) = two * pi * kk(1,i) / rcellx kk(2,i) = two * pi * kk(2,i) / rcelly @@ -3742,9 +3742,9 @@ subroutine readDiagInfo ! BEGIN %%%% ASE printing %%%% ! if (inode==ionode .and. write_ase ) then - inquire(io_ase, opened=test_ase) + inquire(io_ase, opened=test_ase) if ( .not. test_ase ) open(io_ase,file=ase_file, status='old', action='write',& - iostat=stat, position='append') + iostat=stat, position='append') write(io_ase,*) write(io_ase,10) nkp do i=1,nkp @@ -3756,7 +3756,7 @@ subroutine readDiagInfo end if ! ! END %%%% ASE printing %%%% - ! + ! ! Write out k-points if(inode==ionode.AND.iprint_init>2) then write(io_lun,51) nkp @@ -3770,7 +3770,7 @@ subroutine readDiagInfo write (io_lun,'(10x,"Value of kT used for smearing: ",f10.6, a2)') & en_conv * kT, en_units(energy_units) - !****lat<$ + !****lat<$ call stop_backtrace(t=backtrace_timer,who='readDiagInfo') !****lat>$ @@ -3793,9 +3793,9 @@ end subroutine readDiagInfo ! ------------------------------------------------------------------------------ ! Subroutine read_input_aux ! ------------------------------------------------------------------------------ - + !!****f* initial_read/read_input_aux * - !! NAME + !! NAME !! read_input_aux !! USAGE !! call read_input_aux(aux) diff --git a/src/utilities/frame.py b/src/utilities/frame.py index 058265b7d..f0b6d7851 100644 --- a/src/utilities/frame.py +++ b/src/utilities/frame.py @@ -1,12 +1,10 @@ import re -import scipy as sp -from pdb import set_trace +import numpy as np # Regular expressions cell_re = re.compile('cell_vectors(.*?)end cell_vectors', re.M | re.S) stress_re = re.compile('stress_tensor(.*?)end stress_tensor', re.M | re.S) position_re = re.compile('positions(.*?)end positions', re.M | re.S) -position_re = re.compile('positions(.*?)end positions', re.M | re.S) velocity_re = re.compile('velocities(.*?)end velocities', re.M | re.S) force_re = re.compile('forces(.*?)end forces', re.M | re.S) @@ -16,12 +14,12 @@ class Frame: def __init__(self, nat, step): self.step = step self.nat = nat - self.species = sp.zeros(nat) - self.r = sp.zeros((nat, 3), dtype='float') - self.v = sp.zeros((nat, 3), dtype='float') - self.f = sp.zeros((nat, 3), dtype='float') - self.lat = sp.zeros((3, 3), dtype='float') - self.stress = sp.zeros((3, 3), dtype='float') + self.species = np.zeros(nat) + self.r = np.zeros((nat, 3), dtype='float') + self.v = np.zeros((nat, 3), dtype='float') + self.f = np.zeros((nat, 3), dtype='float') + self.lat = np.zeros((3, 3), dtype='float') + self.stress = np.zeros((3, 3), dtype='float') self.ke = 0. self.pe = 0. self.E = 0. @@ -74,4 +72,4 @@ def parse_frame(self, buf): bits.pop(0) bits.pop(0) for j in range(3): - self.f[i,j] = float(bits[j]) \ No newline at end of file + self.f[i,j] = float(bits[j]) diff --git a/src/utilities/md_analysis.py b/src/utilities/md_analysis.py index 636d96cb9..9162f1b20 100755 --- a/src/utilities/md_analysis.py +++ b/src/utilities/md_analysis.py @@ -4,14 +4,14 @@ import sys import re import os.path -import scipy as sp +import numpy as np import matplotlib.pyplot as plt from frame import Frame from md_tools import Pairdist, MSER, VACF, MSD, autocorr -from pdb import set_trace ha2ev = 27.211399 ha2k = 3.15737513e5 +ha_bohr3_to_gpa = 29421.01549104606 # Regular expressions frame_re = re.compile('frame') @@ -61,7 +61,7 @@ def parse_init_config(conf_filename): a = [float(bit) for bit in infile.readline().strip().split()] b = [float(bit) for bit in infile.readline().strip().split()] c = [float(bit) for bit in infile.readline().strip().split()] - data['latvec'] = sp.array([a,b,c]) + data['latvec'] = np.array([a,b,c]) natoms = int(infile.readline().strip()) data['natoms'] = natoms coords = [] @@ -70,8 +70,8 @@ def parse_init_config(conf_filename): x, y, z, spec, cx, cy, cz = infile.readline().strip().split() coords.append([float(x), float(y), float(z)]) species.append(int(spec)) - data['coords'] = sp.array(coords) - data['species'] = sp.array(species) + data['coords'] = np.array(coords) + data['species'] = np.array(species) scount = {} for i in range(natoms): if data['species'][i] in scount.keys(): @@ -84,14 +84,10 @@ def parse_init_config(conf_filename): return data def read_stats(stats_file, nstop): - nstep = 0 data = {} header = True with open(stats_file, 'r') as statfile: for line in statfile: - if nstop != -1: - if nstep > nstop: - break if header: col_id = line.strip().split() for col in col_id: @@ -99,16 +95,18 @@ def read_stats(stats_file, nstop): header = False else: bits = line.strip().split() + step = int(bits[0]) + if nstop != -1 and step > nstop: + break for i, bit in enumerate(bits): if i==0: - info = int(bit) + info = step else: info = float(bit) data[col_id[i]].append(info) - nstep += 1 for key in data: - data[key] = sp.array(data[key]) - return nstep, data + data[key] = np.array(data[key]) + return len(data['step']), data # Command line arguments parser = argparse.ArgumentParser(description='Analyse a Conquest MD \ @@ -118,13 +116,15 @@ def read_stats(stats_file, nstop): in directories specified by -d') parser.add_argument('-d', '--dirs', nargs='+', default=['.',], dest='dirs', action='store', help='Directories to compare') -parser.add_argument('--description', nargs='+', default='', dest='desc', +parser.add_argument('--description', nargs='+', default=[], dest='desc', action='store', help='Description of graph for legend \ (only if using --compare)') parser.add_argument('-f', '--frames', action='store', dest='framesfile', - default='Frames', help='MD frames file') -parser.add_argument('-s', '--stats', action='store', dest='statfile', - default='Stats', help='MD statistics file') + default='md.frames', help='MD frames file') +parser.add_argument('-s', '--stats-file', action='store', dest='statfile', + default='md.stats', help='MD statistics file') +parser.add_argument('--heatflux-file', action='store', dest='heatfluxfile', + default='md.heatflux', help='MD heat-flux file') parser.add_argument('--skip', action='store', dest='nskip', default=0, type=int, help='Number of equilibration steps to skip') parser.add_argument('--stride', action='store', dest='stride', default=1, @@ -153,11 +153,14 @@ def read_stats(stats_file, nstop): help='Generate plot with landscape orientation') parser.add_argument('--pub', action='store_true', dest='pub', help='Publication text size') -parser.add_argument('--nbins', action='store', dest='nbins', default=100, - help='Number of histogram bins') -parser.add_argument('--rdfwidth', action='store', dest='rdfwidth', - default=0.05, help='RDF histogram bin width (A)') +rdf_bins = parser.add_mutually_exclusive_group() +rdf_bins.add_argument('--nbins', action='store', dest='nbins', default=100, + type=int, help='Number of RDF histogram bins') +rdf_bins.add_argument('--rdfwidth', action='store', dest='rdfwidth', + default=None, type=float, + help='RDF histogram bin width (A)') parser.add_argument('--rdfcut', action='store', dest='rdfcut', default=10.0, + type=float, help='Distance cutoff for RDF') parser.add_argument('--dump', action='store_true', dest='dump', help='Dump secondary data used to generate plots') @@ -165,6 +168,9 @@ def read_stats(stats_file, nstop): type=str, help='Compute MSER for the given property') opts = parser.parse_args() +if opts.pub: + plt.rcParams.update({'font.size': 14, 'axes.labelsize': 14, + 'legend.fontsize': 12}) if (opts.vacf or opts.msd or opts.stress or opts.rdf): read_frames = True else: @@ -183,119 +189,116 @@ def read_stats(stats_file, nstop): natoms = init_config['natoms'] dt = float(cq_params['AtomMove.Timestep']) species = cq_params['species'] - extended_system = False - if 'MD.Thermostat' in cq_params.keys(): - if cq_params['MD.Thermostat'] == 'nhc': - extended_system = True - if cq_params['MD.Thermostat'] == 'ssm': - extended_system = True - if cq_params['MD.Thermostat'] == 'svr': - extended_system = False - if 'MD.Barostat' in cq_params.keys(): - if cq_params['MD.Barostat'] == 'iso-ssm': - extended_system = True - if cq_params['MD.Barostat'] == 'ortho-ssm': - extended_system = True - if cq_params['MD.Barostat'] == 'iso-mttk': - extended_system = True - - # Parse the statistics file - nsteps, data = read_stats(opts.statfile,opts.nstop) - avg = {} - std = {} - for key in data: - data[key] = sp.array(data[key]) - avg[key] = sp.mean(data[key][opts.nequil:-1]) - std[key] = sp.std(data[key][opts.nequil:-1]) - time = [float(s)*dt for s in data['step']] - data['time'] = sp.array(time) - - # Plot the statistics - if opts.landscape: - fig1, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2, ncols=2, sharex=True, figsize=(11,7)) - plt.tight_layout(pad=6.5) - else: - fig1, (ax1, ax2, ax3, ax4) = plt.subplots(nrows=4, ncols=1, sharex=True, figsize=(7,10)) - - ax1.plot(data['time'][opts.nskip:], data['pe'][opts.nskip:], 'r-', label='Potential energy') - ax1a = ax1.twinx() - ax1a.plot(data['time'][opts.nskip:], data['ke'][opts.nskip:], 'b-', label='Kinetic energy') - if cq_params['MD.Ensemble'][2] == 't': - if cq_params['MD.Thermostat'] == 'nhc': - ax1a.plot(data['time'][opts.nskip:], data['thermostat'][opts.nskip:], 'g-', label='Thermostat energy') - if cq_params['MD.Thermostat'] == 'svr': - ax1a.plot(data['time'][opts.nskip:], data['thermostat'][opts.nskip:], 'g-', label='Thermostat energy') - if cq_params['MD.Ensemble'][1] == 'p': - if extended_system: - ax1a.plot(data['time'][opts.nskip:], data['box'][opts.nskip:], 'c-', label='Barostat energy') - ax1a.plot(data['time'][opts.nskip:], data['pV'][opts.nskip:], 'm-', label='pV') - ax2.plot(data['time'][opts.nskip:], data['H\''][opts.nskip:]) - ax2.plot((opts.nskip,data['time'][-1]), (avg['H\''],avg['H\'']), '-', - label=r'$\langle H\' \rangle$ = {0:>12.4f} $\pm$ {1:<12.4f}'.format(avg['H\''], std['H\''])) - ax3.plot(data['time'][opts.nskip:], data['T'][opts.nskip:]) - ax3.plot((opts.nskip,data['time'][-1]), (avg['T'],avg['T']), '-', - label=r'$\langle T \rangle$ = {0:>12.4f} $\pm$ {1:<12.4f}'.format(avg['T'], std['T'])) - ax4.plot(data['time'][opts.nskip:], data['P'][opts.nskip:], 'b-') - ax4.plot((opts.nskip,data['time'][-1]), (avg['P'],avg['P']), 'b--', - label=r'$\langle P \rangle$ = {0:>12.4f} $\pm$ {1:<12.4f}'.format(avg['P'], std['P'])) - if cq_params['MD.Ensemble'][1] == 'p': - ax4a = ax4.twinx() - ax4a.plot(data['time'][opts.nskip:], data['V'][opts.nskip:], 'r-') - ax4a.plot((opts.nskip,data['time'][-1]), (avg['V'],avg['V']), 'r--', - label=r'$\langle V \rangle$ = {0:>12.4f} $\pm$ {1:<12.4f}'.format(avg['V'], std['V'])) - ax1.set_ylabel("E (Ha)") - ax2.set_ylabel("H$'$ (Ha)") - ax3.set_ylabel("T (K)") - ax4.set_ylabel("P (GPa)", color='b') - if cq_params['MD.Ensemble'][1] == 'p': - ax4a.set_ylabel("V ($a_0^3$)", color='r') - ax4.set_xlabel("time (fs)") - ax1.legend(loc="upper left") - ax1a.legend(loc="lower right") - ax2.legend() - ax3.legend() - ax4.legend(loc="upper left") - if cq_params['MD.Ensemble'][1] == 'p': - ax4a.legend(loc="lower right") - plt.xlim((opts.nskip,data['time'][-1])) - fig1.subplots_adjust(hspace=0) - fig1.savefig("stats.pdf", bbox_inches='tight') + other_analysis = (opts.vacf or opts.hfacf or opts.msd or opts.rdf or + opts.stress or opts.mser_var) + plot_statistics = opts.stats or not other_analysis + needs_statistics = plot_statistics or opts.mser_var + if needs_statistics: + # Parse the statistics file + nsteps, data = read_stats(opts.statfile,opts.nstop) + avg = {} + std = {} + for key in data: + data[key] = np.array(data[key]) + time = [float(s)*dt for s in data['step']] + data['time'] = np.array(time) + plot_mask = data['step'] >= opts.nskip + equil_mask = data['step'] >= opts.nequil + for key in data: + avg[key] = np.mean(data[key][equil_mask]) + std[key] = np.std(data[key][equil_mask]) + plot_start = data['time'][plot_mask][0] + + if plot_statistics: + # Plot the statistics + if opts.landscape: + fig1, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2, ncols=2, sharex=True, figsize=(11,7)) + plt.tight_layout(pad=6.5) + else: + fig1, (ax1, ax2, ax3, ax4) = plt.subplots(nrows=4, ncols=1, sharex=True, figsize=(7,10)) + + ax1.plot(data['time'][plot_mask], data['pe'][plot_mask], 'r-', label='Potential energy') + ax1a = ax1.twinx() + ax1a.plot(data['time'][plot_mask], data['ke'][plot_mask], 'b-', label='Kinetic energy') + if cq_params['MD.Ensemble'][2] == 't': + if cq_params['MD.Thermostat'] == 'nhc': + ax1a.plot(data['time'][plot_mask], data['thermostat'][plot_mask], 'g-', label='Thermostat energy') + if cq_params['MD.Thermostat'] == 'svr': + ax1a.plot(data['time'][plot_mask], data['thermostat'][plot_mask], 'g-', label='Thermostat energy') + if cq_params['MD.Ensemble'][1] == 'p': + if 'barostat' in data: + ax1a.plot(data['time'][plot_mask], data['barostat'][plot_mask], 'c-', label='Barostat energy') + ax1a.plot(data['time'][plot_mask], data['pV'][plot_mask], 'm-', label='pV') + ax2.plot(data['time'][plot_mask], data['H\''][plot_mask]) + ax2.plot((plot_start,data['time'][-1]), (avg['H\''],avg['H\'']), '-', + label=r'$\langle H\' \rangle$ = {0:>12.4f} $\pm$ {1:<12.4f}'.format(avg['H\''], std['H\''])) + ax3.plot(data['time'][plot_mask], data['T'][plot_mask]) + ax3.plot((plot_start,data['time'][-1]), (avg['T'],avg['T']), '-', + label=r'$\langle T \rangle$ = {0:>12.4f} $\pm$ {1:<12.4f}'.format(avg['T'], std['T'])) + ax4.plot(data['time'][plot_mask], data['P'][plot_mask], 'b-') + ax4.plot((plot_start,data['time'][-1]), (avg['P'],avg['P']), 'b--', + label=r'$\langle P \rangle$ = {0:>12.4f} $\pm$ {1:<12.4f}'.format(avg['P'], std['P'])) + if cq_params['MD.Ensemble'][1] == 'p': + ax4a = ax4.twinx() + ax4a.plot(data['time'][plot_mask], data['V'][plot_mask], 'r-') + ax4a.plot((plot_start,data['time'][-1]), (avg['V'],avg['V']), 'r--', + label=r'$\langle V \rangle$ = {0:>12.4f} $\pm$ {1:<12.4f}'.format(avg['V'], std['V'])) + ax1.set_ylabel("E (Ha)") + ax2.set_ylabel("H$'$ (Ha)") + ax3.set_ylabel("T (K)") + ax4.set_ylabel("P (GPa)", color='b') + if cq_params['MD.Ensemble'][1] == 'p': + ax4a.set_ylabel("V ($a_0^3$)", color='r') + ax4.set_xlabel("time (fs)") + ax1.legend(loc="upper left") + ax1a.legend(loc="lower right") + ax2.legend() + ax3.legend() + ax4.legend(loc="upper left") + if cq_params['MD.Ensemble'][1] == 'p': + ax4a.legend(loc="lower right") + plt.xlim((plot_start,data['time'][-1])) + fig1.subplots_adjust(hspace=0) + fig1.savefig("stats.pdf", bbox_inches='tight') else: # If we're comparing statistics in several directories, use a simplified plot fig1, (ax1, ax2, ax3) = plt.subplots(nrows=3, ncols=1, sharex=True, figsize=(7,7)) ax1a = ax1.twinx() + labels = opts.desc if opts.desc else opts.dirs + if len(labels) != len(opts.dirs): + parser.error('--description must provide one label per comparison directory') + time_limits = [] for ind, d in enumerate(opts.dirs): path = os.path.join(d, cq_input_file) cq_params = parse_cq_input(path) - path = os.path.join(d, cq_params['IO.Coordinates']) - init_config = parse_init_config(path) - natoms = init_config['natoms'] dt = float(cq_params['AtomMove.Timestep']) - species = cq_params['species'] path = os.path.join(d, opts.statfile) nsteps, data = read_stats(path,opts.nstop) time = [float(s)*dt for s in data['step']] - data['time'] = sp.array(time) - - ax1.plot(data['time'][opts.nskip:], data['H\''][opts.nskip:], - linewidth=0.5, label=opts.desc[ind]) - y1,y2 = ax1.get_ylim() - ax1a.set_ylim(y1*ha2k,y2*ha2k) - ax2.plot(data['time'][opts.nskip:], data['T'][opts.nskip:], - linewidth=0.5, label=opts.desc[ind]) - ax3.plot(data['time'][opts.nskip:], data['P'][opts.nskip:], - linewidth=0.5, label=opts.desc[ind]) - - ax1.set_ylabel("H$'$ (Ha)") - ax1a.set_ylabel("H$'$ (K)") - ax2.set_ylabel("T (K)") - ax3.set_ylabel("P (GPa)") - ax3.set_xlabel("time (fs)") - ax1.legend() - plt.xlim((opts.nskip,data['time'][-1])) - fig1.subplots_adjust(hspace=0) - fig1.savefig("stats.pdf", bbox_inches='tight') + data['time'] = np.array(time) + plot_mask = data['step'] >= opts.nskip + time_limits.append((data['time'][plot_mask][0], data['time'][-1])) + + ax1.plot(data['time'][plot_mask], data['H\''][plot_mask], + linewidth=0.5, label=labels[ind]) + ax2.plot(data['time'][plot_mask], data['T'][plot_mask], + linewidth=0.5, label=labels[ind]) + ax3.plot(data['time'][plot_mask], data['P'][plot_mask], + linewidth=0.5, label=labels[ind]) + + y1,y2 = ax1.get_ylim() + ax1a.set_ylim(y1*ha2k,y2*ha2k) + ax1.set_ylabel("H$'$ (Ha)") + ax1a.set_ylabel("H$'$ (K)") + ax2.set_ylabel("T (K)") + ax3.set_ylabel("P (GPa)") + ax3.set_xlabel("time (fs)") + ax1.legend() + plt.xlim((min(limit[0] for limit in time_limits), + max(limit[1] for limit in time_limits))) + fig1.subplots_adjust(hspace=0) + fig1.savefig("stats.pdf", bbox_inches='tight') # Plot MSER if opts.mser_var: @@ -307,38 +310,51 @@ def read_stats(stats_file, nstop): # Plot heat flux autocorrelation function if opts.hfacf: - window = int(opts.acfwindow // dt) - G = sp.zeros((3,3,window)) - time = sp.array([float(i)*dt for i in range(window)]) - nruns = 0 - for ind, d in enumerate(opts.dirs): - path = os.path.join(d, heatfluxfile) - + flux_runs = [] + flux_dt = None + for d in opts.dirs: + run_params = parse_cq_input(os.path.join(d, cq_input_file)) + run_dt = float(run_params['AtomMove.Timestep']) + if flux_dt is None: + flux_dt = run_dt + elif not np.isclose(run_dt, flux_dt): + parser.error('all heat-flux trajectories must use the same timestep') + path = os.path.join(d, opts.heatfluxfile) J = [] - t = [] - nsteps = 0 with open(path, 'r') as infile: for line in infile: step, Jx, Jy, Jz = line.split() step = int(step) - Jx = float(Jx) - Jy = float(Jy) - Jz = float(Jz) - J.append([Jx, Jy, Jz]) - t.append(step*dt) - nsteps += 1 - J = sp.array(J) - t = sp.array(t) - - nwindows = int((nsteps - opts.nskip) // window) - for i in range(3): - for j in range(3): - for k in range(nwindows): - nruns += 1 - start = opts.nskip + k*window - finish = opts.nskip + k*window + window - G += autocorr(J[start:finish,i],J[start:finish,j]) - G = G / float(nruns) + if step < opts.nskip: + continue + if opts.nstop != -1 and step > opts.nstop: + break + J.append([float(Jx), float(Jy), float(Jz)]) + if J: + flux_runs.append(np.array(J)) + + if not flux_runs: + parser.error('no heat-flux samples remain after step selection') + if opts.acfwindow > 0.0: + window = int(opts.acfwindow // flux_dt) + else: + window = min(len(run) for run in flux_runs) + if window < 1 or any(len(run) < window for run in flux_runs): + parser.error('autocorrelation window exceeds the available heat-flux data') + + G = np.zeros((3,3,window)) + nruns = 0 + for J in flux_runs: + for start in range(0, len(J)-window+1, window): + sample = J[start:start+window] + for i in range(3): + for j in range(3): + G[i,j,:] += autocorr(sample[:,i], sample[:,j]) + nruns += 1 + if nruns == 0: + parser.error('heat-flux data do not contain a complete correlation window') + G /= float(nruns) + time = np.arange(window, dtype=float)*flux_dt plt.figure("HFACF") plt.xlabel("t (fs)") plt.ylabel("HFACF") @@ -377,23 +393,26 @@ def read_stats(stats_file, nstop): continue else: done = True - if n <= opts.nskip: + if n < opts.nskip: continue elif n%opts.stride != 0: continue - else: - nframes += 1 - if opts.nstop != -1: - if n > opts.nstop: - done = True + if opts.nstop != -1 and n > opts.nstop: + done = True + continue + nframes += 1 sys.stdout.write("Processing frame {}\r".format(n)) if first_frame: first_frame = False f1 = Frame(natoms,n) f1.parse_frame(buf) if opts.rdf: + if opts.rdfwidth is None: + rdfwidth = opts.rdfcut/opts.nbins + else: + rdfwidth = opts.rdfwidth pairdist = Pairdist(natoms, init_config['nspecies'], - float(opts.rdfcut), float(opts.rdfwidth), + opts.rdfcut, rdfwidth, cq_params['species'], init_config['species_count']) if opts.vacf: @@ -406,8 +425,8 @@ def read_stats(stats_file, nstop): time.append(n*dt) if opts.stress: - stress.append(f.stress) lat.append(f.lat) + stress.append(f.stress/(f.lat[0,0]*f.lat[1,1]*f.lat[2,2])) if opts.rdf: pairdist.update_rdf(f) if opts.vacf: @@ -421,35 +440,29 @@ def read_stats(stats_file, nstop): else: buf += line - time = data['time'] + time = np.array(time) time = time - time[0] print() print("Analysing {} frames...".format(nframes)) # Plot the stress if opts.stress: - stress = sp.array(stress) - lat = sp.array(lat) - mean_stress = sp.zeros((3,3)) - mean_lat = sp.zeros((3,3)) - for i in range(3): - for j in range(3): - mean_stress[i,j] = sp.mean(stress[:,i,j]) - mean_lat[i,j] = sp.mean(lat[:,i,j]) - plt.figure("Stress") - - if cq_params['MD.Ensemble'][1] == "p": + stress = np.array(stress)*ha_bohr3_to_gpa + lat = np.array(lat) + mean_stress = np.mean(stress, axis=0) + mean_lat = np.mean(lat, axis=0) + variable_cell = cq_params['MD.Ensemble'][1] == "p" + + if variable_cell: fig2, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, sharex=True) else: - fig2, (ax1,) = plt.subplots(nrows=1, ncols=1) + fig2, ax1 = plt.subplots(nrows=1, ncols=1) - plt.xlabel("t (fs)") ax1.set_ylabel("Stress (GPa)") - ax2.set_ylabel("Cell dimension ($a_0$)") - plt.xlim((time[opts.nskip], time[-1])) - ax1.plot(time[opts.nskip:], stress[:,0,0], 'r-', label='xx', linewidth=1.0) - ax1.plot(time[opts.nskip:], stress[:,1,1], 'g-', label='yy', linewidth=1.0) - ax1.plot(time[opts.nskip:], stress[:,2,2], 'b-', label='zz', linewidth=1.0) + ax1.set_xlim((time[0], time[-1])) + ax1.plot(time, stress[:,0,0], 'r-', label='xx', linewidth=1.0) + ax1.plot(time, stress[:,1,1], 'g-', label='yy', linewidth=1.0) + ax1.plot(time, stress[:,2,2], 'b-', label='zz', linewidth=1.0) ax1.plot((time[0],time[-1]), (mean_stress[0,0], mean_stress[0,0]), 'r-', label=r'$\langle S_{{xx}} \rangle$ = {0:<10.4f}'.format(mean_stress[0,0])) ax1.plot((time[0],time[-1]), (mean_stress[1,1], mean_stress[1,1]), 'g-', @@ -457,21 +470,24 @@ def read_stats(stats_file, nstop): ax1.plot((time[0],time[-1]), (mean_stress[2,2], mean_stress[2,2]), 'b-', label=r'$\langle S_{{zz}} \rangle$ = {0:<10.4f}'.format(mean_stress[2,2])) - if cq_params['MD.Ensemble'][1] == "p": - ax2.plot(time[opts.nskip:], lat[:,0,0], 'r-', label='a', linewidth=1.0) - ax2.plot(time[opts.nskip:], lat[:,1,1], 'g-', label='b', linewidth=1.0) - ax2.plot(time[opts.nskip:], lat[:,2,2], 'b-', label='c', linewidth=1.0) + ax1.legend(bbox_to_anchor=(1.05,1), loc=2, borderaxespad=0.) + if variable_cell: + ax2.set_ylabel("Cell dimension ($a_0$)") + ax2.set_xlabel("t (fs)") + ax2.plot(time, lat[:,0,0], 'r-', label='a', linewidth=1.0) + ax2.plot(time, lat[:,1,1], 'g-', label='b', linewidth=1.0) + ax2.plot(time, lat[:,2,2], 'b-', label='c', linewidth=1.0) ax2.plot((time[0],time[-1]), (mean_lat[0,0], mean_lat[0,0]), 'r-', label=r'$\langle a \rangle$ = {0:<10.4f}'.format(mean_lat[0,0])) ax2.plot((time[0],time[-1]), (mean_lat[1,1], mean_lat[1,1]), 'g-', label=r'$\langle b \rangle$ = {0:<10.4f}'.format(mean_lat[1,1])) ax2.plot((time[0],time[-1]), (mean_lat[2,2], mean_lat[2,2]), 'b-', label=r'$\langle c \rangle$ = {0:<10.4f}'.format(mean_lat[2,2])) - ax1.legend(bbox_to_anchor=(1.05,1), loc=2, borderaxespad=0.) ax2.legend(bbox_to_anchor=(1.05,1), loc=2, borderaxespad=0.) fig2.subplots_adjust(hspace=0) - plt.setp([a.get_xticklabels() for a in fig1.axes[:-1]], visible=False) - fig2.savefig("stress.pdf", bbox_inches='tight') + else: + ax1.set_xlabel("t (fs)") + fig2.savefig("stress.pdf", bbox_inches='tight') # Plot the rdf if opts.rdf: diff --git a/src/utilities/md_tools.py b/src/utilities/md_tools.py index d59a232e4..2117267a8 100644 --- a/src/utilities/md_tools.py +++ b/src/utilities/md_tools.py @@ -1,25 +1,22 @@ #!/usr/local/bin/python3 -import scipy as sp +import numpy as np import matplotlib.pyplot as plt -from scipy.linalg import norm -from scipy.integrate import cumtrapz +from scipy.integrate import cumulative_trapezoid from scipy.signal import correlate -from scipy import histogram from math import ceil, pi -from frame import Frame -from pdb import set_trace bohr2ang = 0.529177249 small = 1.0e-3 def autocorr(x, y=None): - """Autocorrelation function""" - if y.any(): - result = correlate(x, y, mode='full') - else: - result = correlate(x, x, mode='full') - return result[result.size // 2:] + """Return an unbiased positive-lag auto- or cross-correlation.""" + if y is None: + y = x + if len(x) != len(y): + raise ValueError('correlation inputs must have the same length') + result = correlate(x, y, mode='full')[len(x)-1:] + return result/np.arange(len(x), 0, -1) def diff_mic(pos1, pos2, cell): """Minimum image convention relative vector (orthorhombic cell only)""" @@ -44,73 +41,86 @@ def __init__(self, nat, nspec, rcut, binwidth, species, species_count): self.nspec = nspec self.rcut = rcut self.binwidth = binwidth - self.nbins = ceil(rcut/binwidth)+1 + self.nbins = ceil(rcut/binwidth) self.spec_count = species_count self.species = species - self.bins = [] - for i in range(self.nbins): - self.bins.append((float(i)*binwidth + binwidth/2.)) - self.bins = sp.array(self.bins) - self.dt = sp.zeros((self.nat,self.nat), dtype='float') - self.freq_total = sp.zeros(self.nbins, dtype='int') - self.freq = sp.zeros((self.nbins,self.nspec,self.nspec), dtype='int') - self.nfac_total = sp.zeros(self.nbins, dtype='float') - self.nfac = sp.zeros((self.nbins,self.nspec,self.nspec), dtype='float') - self.gr_total = sp.zeros(self.nbins, dtype='float') - self.gr = sp.zeros((self.nbins,self.nspec,self.nspec), dtype='float') + self.edges = np.minimum(np.arange(self.nbins+1)*binwidth, rcut) + self.edges[-1] = rcut + self.bins = 0.5*(self.edges[:-1] + self.edges[1:]) + self.inv_volume_sum = 0.0 + self.dt = np.zeros((self.nat,self.nat), dtype='float') + self.freq_total = np.zeros(self.nbins, dtype='int') + self.freq = np.zeros((self.nbins,self.nspec,self.nspec), dtype='int') + self.nfac_total = np.zeros(self.nbins, dtype='float') + self.nfac = np.zeros((self.nbins,self.nspec,self.nspec), dtype='float') + self.gr_total = np.zeros(self.nbins, dtype='float') + self.gr = np.zeros((self.nbins,self.nspec,self.nspec), dtype='float') def update_rdf(self, frame): self.nframes += 1 - cell = sp.zeros(3) + cell = np.zeros(3) for i in range(3): cell[i] = frame.lat[i,i] - self.volume = cell[0]*cell[1]*cell[2]*bohr2ang**3 + self.volume = abs(np.linalg.det(frame.lat))*bohr2ang**3 self.rho = float(self.nat)/self.volume + self.inv_volume_sum += 1.0/self.volume for i in range(self.nat): for j in range(i+1, self.nat): diff = diff_mic(frame.r[i,:], frame.r[j,:], cell)*bohr2ang - self.dt[i,j] = norm(diff) - self.dt[j,i] = norm(diff) + self.dt[i,j] = np.linalg.norm(diff) + self.dt[j,i] = np.linalg.norm(diff) if self.dt[i,j] < self.rcut: - ind = int(round((self.dt[i,j]+self.binwidth)/self.binwidth))-1 + ind = np.searchsorted(self.edges, self.dt[i,j], side='right')-1 self.freq_total[ind] += 2 if self.nspec > 1: - for ispec in range(self.nspec): - for jspec in range(ispec, self.nspec): - if (ispec == frame.species[i]-1 and jspec == frame.species[j]-1): - self.freq[ind, ispec, jspec] += 2 + ispec, jspec = sorted((int(frame.species[i])-1, + int(frame.species[j])-1)) + if ispec == jspec: + self.freq[ind, ispec, jspec] += 2 + else: + self.freq[ind, ispec, jspec] += 1 def norm_rdf(self): """Normalise the RDF""" - const1 = 4.0*pi*(self.binwidth**3)/3.0 - const2 = self.rho*self.nat*self.nframes for i in range(self.nbins): - vshell = (float(i+1)**3 - float(i)**3)*const1 - self.nfac_total[i] = vshell*const2 + vshell = 4.0*pi*(self.edges[i+1]**3-self.edges[i]**3)/3.0 + self.nfac_total[i] = (vshell*self.nat*(self.nat-1)* + self.inv_volume_sum) if self.nspec > 1: for ispec in range(self.nspec): for jspec in range(self.nspec): - const3 = self.rho*self.spec_count[ispec+1]*self.spec_count[jspec+1]/self.nat - self.nfac[i,ispec,jspec] = vshell*const3*self.nframes - self.gr_total = self.freq_total.astype(float)/self.nfac_total + ni = self.spec_count[ispec+1] + nj = self.spec_count[jspec+1] + npairs = ni*(ni-1) if ispec == jspec else ni*nj + self.nfac[i,ispec,jspec] = (vshell*npairs* + self.inv_volume_sum) + self.gr_total = np.divide(self.freq_total, self.nfac_total, + out=np.zeros_like(self.nfac_total), + where=self.nfac_total > 0) if self.nspec > 1: - self.gr = self.freq.astype(float)/self.nfac + self.gr = np.divide(self.freq, self.nfac, + out=np.zeros_like(self.nfac), + where=self.nfac > 0) def get_coordination(self): """Compute coordination""" gxrsq = self.gr_total*self.bins**2 - self.coord_total = sp.zeros(self.nbins, dtype='float') - self.coord_total[1:] = cumtrapz(gxrsq,self.bins) - self.coord_total *= 4.*pi*self.rho + self.coord_total = np.zeros(self.nbins, dtype='float') + self.coord_total[1:] = cumulative_trapezoid(gxrsq,self.bins) + mean_inv_volume = self.inv_volume_sum/self.nframes + self.coord_total *= 4.*pi*(self.nat-1)*mean_inv_volume if self.nspec > 1: - self.coord = sp.zeros((self.nbins,self.nspec,self.nspec), dtype='float') + self.coord = np.zeros((self.nbins,self.nspec,self.nspec), dtype='float') for ispec in range(self.nspec): for jspec in range(ispec,self.nspec): gxrsq = self.gr[:,ispec,jspec]*self.bins**2 - self.coord[1:,ispec,jspec] = cumtrapz(gxrsq[:], self.bins) - self.coord *= 4.*pi*self.rho # check this + self.coord[1:,ispec,jspec] = cumulative_trapezoid(gxrsq[:], self.bins) + nneighbours = self.spec_count[jspec+1] + if ispec == jspec: + nneighbours -= 1 + self.coord[:,ispec,jspec] *= 4.*pi*nneighbours*mean_inv_volume def plot_gr(self): plt.figure("RDF") @@ -120,9 +130,9 @@ def plot_gr(self): else: fig3, axl = plt.subplots() axl.minorticks_on() - axl.grid(b=True, which='major', axis='x', color='gray', linestyle='-') - axl.grid(b=True, which='minor', axis='x', color='gray', linestyle='--') - axl.grid(b=True, which='major', axis='y', color='gray', linestyle='-') + axl.grid(visible=True, which='major', axis='x', color='gray', linestyle='-') + axl.grid(visible=True, which='minor', axis='x', color='gray', linestyle='--') + axl.grid(visible=True, which='major', axis='y', color='gray', linestyle='-') # axl.grid(b=True, which='minor', axis='y', color='gray', linestyle='--') axr = axl.twinx() axl.set_ylabel("g(r)", color='b') @@ -134,9 +144,9 @@ def plot_gr(self): axr.set_ylim(bottom=axl.get_ylim()[0], top=axl.get_ylim()[1]*10.0) if self.nspec > 1: ax2.minorticks_on() - ax2.grid(b=True, which='major', axis='x', color='gray', linestyle='-') - ax2.grid(b=True, which='minor', axis='x', color='gray', linestyle='--') - ax2.grid(b=True, which='major', axis='y', color='gray', linestyle='-') + ax2.grid(visible=True, which='major', axis='x', color='gray', linestyle='-') + ax2.grid(visible=True, which='minor', axis='x', color='gray', linestyle='--') + ax2.grid(visible=True, which='major', axis='y', color='gray', linestyle='-') # ax2.grid(b=True, which='minor', axis='y', color='gray', linestyle='--') for ispec in range(self.nspec): for jspec in range(ispec,self.nspec): @@ -180,12 +190,12 @@ def dump_gr(self): def get_bondlength(self, bondcut, frame, printall): - bond_tot = sp.zeros((self.nspec, self.nspec), dtype=float) - bondsq_tot = sp.zeros((self.nspec, self.nspec), dtype=float) - bond_avg = sp.zeros((self.nspec, self.nspec), dtype=float) - bond_sd = sp.zeros((self.nspec, self.nspec), dtype=float) - bond_min = sp.zeros((self.nspec, self.nspec), dtype=float) - nbonds = sp.zeros((self.nspec, self.nspec), dtype=int) + bond_tot = np.zeros((self.nspec, self.nspec), dtype=float) + bondsq_tot = np.zeros((self.nspec, self.nspec), dtype=float) + bond_avg = np.zeros((self.nspec, self.nspec), dtype=float) + bond_sd = np.zeros((self.nspec, self.nspec), dtype=float) + bond_min = np.zeros((self.nspec, self.nspec), dtype=float) + nbonds = np.zeros((self.nspec, self.nspec), dtype=int) bond_min = bondcut for i in range(self.nat): @@ -209,7 +219,7 @@ def get_bondlength(self, bondcut, frame, printall): for j in range(i,self.nspec): if nbonds[i,j] > 0: bond_avg[i,j] = bond_tot[i,j]/float(nbonds[i,j]) - bond_sd[i,j] = sp.sqrt(bondsq_tot[i,j]/nbonds[i,j] - bond_avg[i,j]**2) + bond_sd[i,j] = np.sqrt(bondsq_tot[i,j]/nbonds[i,j] - bond_avg[i,j]**2) pair = "{}-{}".format(self.species[i+1], self.species[j+1]) print(f'{pair}: {bond_avg[i,j]:>8.4f} +/- {bond_sd[i,j]:>8.4f}') @@ -224,36 +234,36 @@ class MSER: --- K P White, Simulation 69, 323 (1997)""" def __init__(self, nframes, varname, var_traj): - self.n_j = nframes-1 self.propname = varname - self.traj = var_traj - self.mser = sp.zeros(self.n_j, dtype='float') + self.traj = np.asarray(var_traj) + self.n_j = len(self.traj) + if self.n_j < 2: + raise ValueError('MSER requires at least two trajectory samples') + self.mser = np.zeros(self.n_j, dtype='float') # stop before the end otherwise the MSER becomes very noisy - self.mser_cut = 200 + self.mser_cut = min(200, max(1, self.n_j//2)) + self.n_candidates = self.n_j-self.mser_cut def get_point(self, d_j): - prefac = 1.0/(self.n_j-d_j)**2 - ybar_ij = sp.mean(self.traj[d_j:]) - variance = 0.0 - for i in range(d_j+1,self.n_j): - variance += (self.traj[i] - ybar_ij)**2 - return prefac*variance + sample = self.traj[d_j:] + return np.sum((sample-np.mean(sample))**2)/len(sample)**2 def get_mser(self): for i in range(self.n_j): self.mser[i] = self.get_point(i) def mser_min(self): - return sp.argmin(self.mser[:-self.mser_cut]) + return np.argmin(self.mser[:self.n_candidates]) def plot_mser(self, steps): plt.figure("{} MSER".format(self.propname)) plt.xlabel("step") plt.ylabel("MSER ({})".format(self.propname)) - plt.plot(steps[:-200], self.mser[:-200], 'k-') + plt.plot(steps[:self.n_candidates], self.mser[:self.n_candidates], 'k-') mser_min = self.mser_min() - lab = "Minimum at step {}".format(mser_min) - plt.axvline(x=mser_min, label=lab) + minimum_step = steps[mser_min] + lab = "Minimum at step {}".format(minimum_step) + plt.axvline(x=minimum_step, label=lab) plt.legend(loc="upper left") plt.savefig("mser.pdf", bbox_inches='tight') @@ -280,11 +290,11 @@ def update_vacf(self, step, frame): self.steps.append(step) self.vacf.append(0.0) for i in range(self.nat): - self.vacf[-1] += sp.dot(self.init_v[i,:], frame.v[i,:]) + self.vacf[-1] += np.dot(self.init_v[i,:], frame.v[i,:]) def norm_vacf(self): - self.vacf = sp.array(self.vacf)/self.nat - self.time = sp.array(self.steps, dtype='float')*self.dt + self.vacf = np.array(self.vacf)/self.nat + self.time = (np.array(self.steps, dtype='float') - self.steps[0])*self.dt def plot_vacf(self): filename = "vacf.pdf" @@ -311,30 +321,30 @@ def __init__(self, nat, dt, init_frame): self.nat = nat self.dt = dt self.init_r = init_frame.r - self.r_prev = sp.copy(self.init_r) + self.r_prev = np.copy(self.init_r) self.msd = [] self.steps = [] - self.init_cell = sp.zeros(3, dtype='float') - self.r_diff = sp.zeros((self.nat,3), dtype='float') + self.init_cell = np.zeros(3, dtype='float') + self.r_diff = np.zeros((self.nat,3), dtype='float') for i in range(3): self.init_cell[i] = init_frame.lat[i,i] def update_msd(self, step, frame): self.steps.append(step) self.nframes += 1 - cell = sp.zeros(3, dtype='float') + cell = np.zeros(3, dtype='float') for i in range(3): cell[i] = frame.lat[i,i] self.msd.append(0.0) for i in range(self.nat): diff = diff_mic(frame.r[i,:], self.r_prev[i,:], cell) # is this right? self.r_diff[i,:] += diff - self.msd[-1] += sp.sum(self.r_diff[i,:]**2) - self.r_prev = sp.copy(frame.r) + self.msd[-1] += np.sum(self.r_diff[i,:]**2) + self.r_prev = np.copy(frame.r) def norm_msd(self): - self.msd = sp.array(self.msd)/self.nat - self.time = sp.array(self.steps, dtype='float')*self.dt + self.msd = np.array(self.msd)/self.nat + self.time = (np.array(self.steps, dtype='float') - self.steps[0])*self.dt def plot_msd(self): filename = "msd.pdf"