Clarification for the restart/restore or post analysis based on a converged *.nc

[zhangguangping]

Dear all, I have learning on how to restart/resore or even do the post analysis based on a converged *.nc file in ATK-2015.1. I found here the tutorial very useful for the psot analysis based on a converged *.nc file http://www.quantumwise.com/publications/tutorials/item/504-computing-quantities-from-converged-calculations The trick to use a converged *.nc file doing a post analysis with out any futher SCFs is to use nlread, eg., the following is to do a total energy calculation based on the already converged Au-C6H4S2-Au.nc file.

Code

# -------------------------------------------------------------
# Analysis from File
# -------------------------------------------------------------
configuration = nlread('Au-C6H4S2-Au.nc', object_id='gID000')[0]

# -------------------------------------------------------------
# Total Energy
# -------------------------------------------------------------
total_energy = TotalEnergy(configuration)
nlsave('analysis.nc', total_energy)
nlprint(total_energy)

Code

Next, check that the Object id corresponds to the self-consistent calculation you want to base the analysis on. If your calculation was a "standard run", then most likely the default gID001 is the one you want. However, if you did a geometry optimization, gID000 corresponds to the initial geometry, and the optimized geometry will be gID001. You can inspect the contents of the NetCDF file, and the object ids, in the Result Browser in the main VNL job.

Note: nlread() can in principle read many objects at once from a NetCDF file. Object ids must however be unique, so by using the object_id keyword, you are guaranteed to get only one object back. However, nlread() always returns a list of objects, even if there is only one match, so we need the [0] to pick out the first/only configuration matching the object_id.

However, I think there maybe something worth to emphasize. That is,

Code

configuration = nlread('Au-C6H4S2-Au.nc', object_id='gID000')[0]

we should specify the object_id properly and may not be the gID000 one. It depends on your number of nlsave used in the job generating this *.nc file and also the quantities you want to base for the post analysis, usually the object_id contains the converged electron density not others would be expected to used for the post analysis. The number of object_id in the *.nc file depends on the nlsave used in the job generating this *.nc. So, usually and most probably the first object_id that is gID000 is not the one contains the converged electron density. For example, the following *.py input file contains 6 nlsave sentances and the second object_id that is gID001 is what we need to do the total energy post analysis.

Code

# -------------------------------------------------------------
# Two-probe Configuration
# -------------------------------------------------------------

# -------------------------------------------------------------
# Left Electrode
# -------------------------------------------------------------

# Set up lattice
vector_a = [8.65127, 0.0, 0.0]*Angstrom
vector_b = [-4.32564, 7.49222, 0.0]*Angstrom
vector_c = [0.0, 0.0, 7.06373620597]*Angstrom
left_electrode_lattice = UnitCell(vector_a, vector_b, vector_c)

# Define elements
left_electrode_elements = [...
                              ...
                              ...]

# Define coordinates
left_electrode_coordinates = [...
                              ...
                              ...]*Angstrom

# Set up configuration
left_electrode = BulkConfiguration(
    bravais_lattice=left_electrode_lattice,
    elements=left_electrode_elements,
    cartesian_coordinates=left_electrode_coordinates
    )

# -------------------------------------------------------------
# Right Electrode
# -------------------------------------------------------------

# Set up lattice
vector_a = [8.65127, 0.0, 0.0]*Angstrom
vector_b = [-4.32564, 7.49222, 0.0]*Angstrom
vector_c = [0.0, 0.0, 7.06373620597]*Angstrom
right_electrode_lattice = UnitCell(vector_a, vector_b, vector_c)

# Define elements
right_electrode_elements = [...
                              ...
                              ...]

# Define coordinates
right_electrode_coordinates = [...
                              ...
                              ...]*Angstrom

# Set up configuration
right_electrode = BulkConfiguration(
    bravais_lattice=right_electrode_lattice,
    elements=right_electrode_elements,
    cartesian_coordinates=right_electrode_coordinates
    )

# -------------------------------------------------------------
# Central Region
# -------------------------------------------------------------

# Set up lattice
vector_a = [8.65127, 0.0, 0.0]*Angstrom
vector_b = [-4.32564, 7.49222, 0.0]*Angstrom
vector_c = [0.0, 0.0, 26.1041296975]*Angstrom
central_region_lattice = UnitCell(vector_a, vector_b, vector_c)

# Define elements
central_region_elements = [...
                              ...
                              ...]

# Define coordinates
central_region_coordinates = [...
                              ...
                              ...]

# Set up configuration
central_region = BulkConfiguration(
    bravais_lattice=central_region_lattice,
    elements=central_region_elements,
    cartesian_coordinates=central_region_coordinates
    )

device_configuration = DeviceConfiguration(
    central_region,
    [left_electrode, right_electrode]
    )
nlsave('Au-C6H4S2-Au.nc', device_configuration)     # the 1st one  gID000
nlprint(device_configuration)

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Basis Set
#----------------------------------------
basis_set = [
    GGABasis.Hydrogen_DoubleZetaPolarized(element=Hydrogen,filter_mesh_cutoff=200*Rydberg),
    GGABasis.Carbon_DoubleZetaPolarized(element=Carbon,filter_mesh_cutoff=200*Rydberg),
    GGABasis.Sulfur_DoubleZetaPolarized(element=Sulfur,filter_mesh_cutoff=200*Rydberg),
    GGABasis.Gold_SingleZetaPolarized(element=Gold,filter_mesh_cutoff=200*Rydberg),
    ]

#----------------------------------------
# Exchange-Correlation
#----------------------------------------
exchange_correlation = GGA.PBE

#----------------------------------------
# Numerical Accuracy Settings
#----------------------------------------
left_electrode_numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=(6, 6, 20),
    density_mesh_cutoff=200.0*Rydberg,
    )

right_electrode_numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=(6, 6, 20),
    density_mesh_cutoff=200.0*Rydberg,
    )

device_numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=(6, 6, 20),
    density_mesh_cutoff=200.0*Rydberg,
    )

#----------------------------------------
# Iteration Control Settings
#----------------------------------------
left_electrode_iteration_control_parameters = IterationControlParameters(
    max_steps=500,
    )

right_electrode_iteration_control_parameters = IterationControlParameters(
    max_steps=500,
    )

device_iteration_control_parameters = IterationControlParameters(
    max_steps=500,
    )

#----------------------------------------
# Poisson Solver Settings
#----------------------------------------
left_electrode_poisson_solver = FastFourier2DSolver(
    boundary_conditions=[[PeriodicBoundaryCondition,PeriodicBoundaryCondition],
                         [PeriodicBoundaryCondition,PeriodicBoundaryCondition],
                         [PeriodicBoundaryCondition,PeriodicBoundaryCondition]]
    )

right_electrode_poisson_solver = FastFourier2DSolver(
    boundary_conditions=[[PeriodicBoundaryCondition,PeriodicBoundaryCondition],
                         [PeriodicBoundaryCondition,PeriodicBoundaryCondition],
                         [PeriodicBoundaryCondition,PeriodicBoundaryCondition]]
    )

#----------------------------------------
# Contour Integral Settings
#----------------------------------------
contour_parameters = DoubleContourIntegralParameters(
    integral_lower_bound=50.0*eV,
    )

#----------------------------------------
# Electrode Calculators
#----------------------------------------
left_electrode_calculator = LCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=left_electrode_numerical_accuracy_parameters,
    iteration_control_parameters=left_electrode_iteration_control_parameters,
    poisson_solver=left_electrode_poisson_solver,
    )

right_electrode_calculator = LCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=right_electrode_numerical_accuracy_parameters,
    iteration_control_parameters=right_electrode_iteration_control_parameters,
    poisson_solver=right_electrode_poisson_solver,
    )

#----------------------------------------
# Device Calculator
#----------------------------------------
calculator = DeviceLCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=device_numerical_accuracy_parameters,
    iteration_control_parameters=device_iteration_control_parameters,
    contour_parameters=contour_parameters,
    electrode_calculators=
        [left_electrode_calculator, right_electrode_calculator],
    )

device_configuration.setCalculator(calculator)
nlprint(device_configuration)
device_configuration.update()
nlsave('Au-C6H4S2-Au.nc', device_configuration)    # the 2nd one  gID001
nlprint( )

# -------------------------------------------------------------
# Transmission Spectrum
# -------------------------------------------------------------
transmission_spectrum = TransmissionSpectrum(
    configuration=device_configuration,
    energies=numpy.linspace(-2,2,400)*eV,
    kpoints=MonkhorstPackGrid(6,6),
    energy_zero_parameter=AverageFermiLevel,
    infinitesimal=1.36057e-05*eV,
    self_energy_calculator=RecursionSelfEnergy(),
    )
nlsave('Au-C6H4S2-Au.nc', transmission_spectrum)    # the 3rd one gID002
nlprint(transmission_spectrum)

# -------------------------------------------------------------
# Molecular Energy Spectrum
# -------------------------------------------------------------
molecular_energy_spectrum = MolecularEnergySpectrum(
    configuration=device_configuration,
    energy_zero_parameter=FermiLevel,
    projection_list=ProjectionList(elements=[Carbon, Hydrogen, Sulfur])
    )
nlsave('Au-C6H4S2-Au.nc', molecular_energy_spectrum)   # the 4th one gID003
nlprint(molecular_energy_spectrum)

# -------------------------------------------------------------
# Eigenstate
# -------------------------------------------------------------
eigenstate = Eigenstate(
    configuration=device_configuration,
    projection_list=ProjectionList(elements=[Carbon, Hydrogen, Sulfur]),
    quantum_number=19,
    )
nlsave('Au-C6H4S2-Au.nc', eigenstate)                             # the 5th one   gID004
nlprint( )

# -------------------------------------------------------------
# Eigenstate
# -------------------------------------------------------------
eigenstate = Eigenstate(
    configuration=device_configuration,
    projection_list=ProjectionList(elements=[Carbon, Hydrogen, Sulfur]),
    quantum_number=20,
    )
nlsave('Au-C6H4S2-Au.nc', eigenstate)                            # the 6th one   gID005
nlprint( )

The first object_id, i.e., gID000 only contain only the configuration (lattice vectors and atom positions) of the two probe system. The converged electron density is stored in the second object_id, i.e., gID001. So, the post analysis input file should read

Code

# -------------------------------------------------------------
# Analysis from File
# -------------------------------------------------------------
configuration = nlread('Au-C6H4S2-Au.nc', object_id='gID001')[0]

# -------------------------------------------------------------
# Total Energy
# -------------------------------------------------------------
total_energy = TotalEnergy(configuration)
nlsave('analysis.nc', total_energy)
nlprint(total_energy)

[Anders Blom]

It simply depends on whether you read the whole file, or parts of it. The following points should be noted:

nlread always returns a list, even if the file only has one object, or the selection criteria only give one object. So you always need to select the object in the list, like [ 0 ] or [-1] etc.

The indexing [1] or [4] etc is NOT for the file, it's for the list returned by nlread. So yes, if you read the whole file
x=nlread("Au-C6H4S2-Au.nc')
then your molecular eigenspectrum is x[3].

But if you only read one object, like
x=nlread('Au-C6H4S2-Au.nc', object_id="gID003")
then the object, whatever it is, is x[0] (or x[-1], which is the same thing in this case).

You can also read all objects of a particular type, like
x=nlread('Au-C6H4S2-Au.nc', class_type=Eigenstate)
Now, x[ 0 ] is the first eigenstate from the NC file, x[1] the second, and x[-1] the last.

All this is basically just regular Python, nothing special we made, as long as you note that, as mentioned, nlread always returns a list, and the keywords object_id and class_type are used to select what to read from the file.

NOTE: I write [ 0 ] with a space because without it, the Forum renders [0] in a funny way. In the code, you can write [0].

[zhangguangping]

Yah, an excellent explaination of this! Indeed, when I post the thread, I realized that the [ 0 ] is for counting the list. I should mean gD000 or gD001, etc.

Thank you very much Anders.

With best regards,

/Guangping

[Anders Blom]

Also note that you can name your objects as you like, in nlsave

nlsave("file.nc", bulk_configuration, object_id="relaxed config")

This can make it easier to use nlread later. However, if you do give object_id, then if an object with that name already exists in the NC file, it will be replaced (without warning). When you do nlsave without it, the new object is always just appended to the file.

[zhangguangping]

Yeh, good to know this.

Thanks so much.

[zhangguangping]

It simply depends on whether you read the whole file, or parts of it. The following points should be noted:

Dear Anders, For saving time by initializing a calculation from the converged state of another, for example, use the converged 0.0V *nc file to initialized the 1.0V calcualtion for the very same system except the bias increased from 0.0 V to 1.0 V, as described in this tutorial http://www.quantumwise.com/publications/tutorials/item/503-saving-time-by-initializing-a-calculation-from-the-converged-state-of-another if I suppose the core part of the input file for 0.0 V is as follows

Code

#----------------------------------------
# Device Calculator
#----------------------------------------
calculator = DeviceLCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=device_numerical_accuracy_parameters,
    iteration_control_parameters=device_iteration_control_parameters,
    contour_parameters=contour_parameters,
    electrode_calculators=
        [left_electrode_calculator, right_electrode_calculator],
    )

device_configuration.setCalculator(calculator)
nlprint(device_configuration)
device_configuration.update()
nlsave('Au-C6H4S2-Au-0.0.nc', device_configuration)

I need to reivse this part of the input file into

Code

#----------------------------------------
# Device Calculator
#----------------------------------------
calculator = DeviceLCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=device_numerical_accuracy_parameters,
    electrode_calculators=
        [left_electrode_calculator, right_electrode_calculator],
    electrode_voltages=( 0.0*Volt, 1.0*Volt)
    )

device_configuration.setCalculator(calculator)

# -------------------------------------------------------------
# Initial State
# -------------------------------------------------------------
old_calculation = nlread('Au-C6H4S2-Au-0.0.nc', DeviceConfiguration)[0]
device_configuration.setCalculator(
    calculator,
    initial_state=old_calculation,
)
device_configuration.update()
nlsave('Au-C6H4S2-Au-1.0.nc', device_configuration)
nlprint(device_configuration)

or even a more simpler one

Code

# -------------------------------------------------------------
# Initial State
# -------------------------------------------------------------
old_calculation = nlread('Au-C6H4S2-Au-0.0.nc', DeviceConfiguration)[0]
device_configuration.setCalculator(
    calculator(electrode_voltages=( 0.0*Volt, 1.0*Volt)),
    initial_state=old_calculation,
)
device_configuration.update()
nlsave('Au-C6H4S2-Au-1.0.nc', device_configuration)
nlprint(device_configuration)

Am I correct? As I understand, device_configuration.setCalculator() does not trigger a real calculation, but just update the settings for the calculation, while device_configuration.update() indeed would trigger the DFT-NEGF SCFs, right? For the above script for 1.0 V, the exectuable would do the left and right electrode calculations again? If so, how to aviod a repeat calculation for the electrodes here, since only shifts of the chemcialpotential are applied to the electrodes? I do not find a setCalculator method for LCAOCalculator used for electrode calculations as that for DeviceLCAOCalculator, where an initial_state method can be used to give an initial state for the electrode calcualtions. Also, I did not find the I-V mini-tutorial  mentioned in the above tutorial for stepping up in bias in I-V calculations. With respect to electrode_voltages=( 0.0*Volt, 1.0*Volt), how would ATK do for the shift of chemicalpotentials of the electrodes relative to the equilibrium ones? Will ATK only shift the right electrode down by 1.0 eV? Or instead, ATK would shift down the right electrode by 0.5 eV and shift up the left electrode by 0.5 eV. I think the latter technique would be favorable to the convergence. With best regards, /Guangping

[Anders Blom]

You are correct, but no, it will not redo the electrodes just because you change the voltage.

That mini-tutorial became a bit outdated so it's been removed. The basic idea of it was however exactly what you have already shown in your scripts, to loop in bias using the lower bias point as initial_state.

A more "modern" take on the same topic is presented here:
http://docs.quantumwise.com/tutorials/atk_transport_calculations.html

[zhangguangping]

Dear Anders,

Thanks for your reply very much.

With respect to electrode_voltages=( 0.0*Volt, 1.0*Volt), how would ATK do for the shift of chemicalpotentials of the electrodes relative to the equilibrium ones? Will ATK only shift the right electrode down by 1.0 eV? Or instead, ATK would shift down the right electrode by 0.5 eV and shift up the left electrode by 0.5 eV. I think the latter technique would be favorable to the convergence.

(1) So, how about the deal with the voltage applized to the electrodes in ATK?

(2) And I am puzzling about the force_restart=True in device_configuration.update(), that is what is the difference between device_configuration.update() and device_configuration.update(force_restart=True) since they both trigger a calculation?

(3) Can you please give more instructions on what is the idea of giving an initial guess for the configuration with anti-parallel spin polarizaton of the electrodes in a magnetic tunnel junction (MTJ) from the converged state of the parallel electrode calculation.Special example: anti-parallel spin in MTJs?
http://www.quantumwise.com/publications/tutorials/item/503-saving-time-by-initializing-a-calculation-from-the-converged-state-of-another
Because I think the convergence behavior and results of spin polarization calculations depend severly on the initial guess. So, a large deviation of parallel electrode calculation from an anti-parallel spin polarizaton of the electrodes would be good initial guess for the anti-parallel spin polarizaton of the electrodes? I am really hard to image that.

With best regards,

/Guangping

[zhangguangping]

Dear Anders, For the 2nd qestion, I read the tutorial here, and get some idea of the difference between device_configuration.update() and device_configuration.update(force_restart=True). http://www.quantumwise.com/publications/tutorials/item/502-restarting-stopped-calculations As I understand, device_configuration.update() is used with nlread() which give an initial_state and resulting a calculation with a previous result as an initial guess, this can also a restart of a previous stopped calculation. For example:

Code

old_calculation = nlread('Au-C6H4S2-Au-0.0.nc', DeviceConfiguration)[0]
device_configuration.setCalculator(
    calculator,
    initial_state=old_calculation,
)
device_configuration.update()
nlsave('Au-C6H4S2-Au-1.0.nc', device_configuration)
nlprint(device_configuration)

However, device_configuration.update(force_restart=True) is used with nlread() which read all the calcualtion settings of a previous calculation, and using the SCF information usually the last SCF density matrix to do a restart calculation of the stopped job. For example, the above script if use the restart method to do the continue calculation should be

Code

device_configuration = nlread('Au-C6H4S2-Au-0.0-checkpoint.nc')[0]
device_configuration.update(force_restart=True)
nlsave('Au-C6H4S2-Au-1.0.nc', device_configuration)
nlprint(device_configuration)

Am I right, if any misunderstanding, please correct them ! Can Au-C6H4S2-Au-0.0.nc which generated by nlsave('Au-C6H4S2-Au-0.0.nc', device_configuration) be used as the checkpoint file for Au-C6H4S2-Au-0.0-checkpoint.nc which should be generated by CheckpointHandler? I think they should be different. The former only contains the converged scf density matrix if the NEGF is done, and if the exectuable stops at the SCF iteration, then there is no density matrix in Au-C6H4S2-Au-0.0.nc. Intead, the latter one should be generated by the CheckpointHandler in DeviceLCAOCalculator for a two probe system for example, which would update itself during the self-consistent calculation, and contains the density matrix of the last scf step. Therefore, Au-C6H4S2-Au-0.0-checkpoint.nc can be used as a restart checkpoint and also as an initial guess for a calculation. If the NEGF is done, Au-C6H4S2-Au-0.0.nc can be used in the same way as Au-C6H4S2-Au-0.0-checkpoint.nc, while if the executable stops during the SCF interation, Au-C6H4S2-Au-0.0.nc either can not used as a restart checkpoint nor as an initial guess for a calculation since there is no density matrix information. Thanks so much ! With best. /Guangping

[Anders Blom]

You are exactly correct on update w/out "force".

You should absolutely use a separate checkpoint file compared to whatever you nlsave into! And most of the time, the checkpoint file is not a good starting guess for another calculation since it will not be the last step, only the last checkpoint step which may be step 9 if you converged at step 11.

[zhangguangping]

Dear Anders, Thanks very much for your kind reply. For example the following script,

Code

old_calculation = nlread('Au-C6H4S2-Au-0.0.nc', DeviceConfiguration)[0]
device_configuration.setCalculator(
    calculator,
    initial_state=old_calculation,
)
device_configuration.update()
nlsave('Au-C6H4S2-Au-1.0.nc', device_configuration)
nlprint(device_configuration)

the Au-C6H4S2-Au-0.0.nc file read by nlread must be a checkpoint file generated by CheckpointHandler not written by nlsave? And the sentance nlsave('Au-C6H4S2-Au-1.0.nc', device_configuration) above does not write the converged density matrix into Au-C6H4S2-Au-0.0.nc even if it follows the self-consistent scf calculation device_configuration.update() and can not used as an initial guess for another calculation ? As I tested, it is not the above case, that is  Au-C6H4S2-Au-0.0.nc written by nlsave can be used as a checkpoint file to at least restart the same calculation. So in this case, it can also used as an initial guess for another calcualtion. The test script is,

Code

#----------------------------------------
# Device Calculator
#----------------------------------------
calculator = DeviceLCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=device_numerical_accuracy_parameters,
    iteration_control_parameters=device_iteration_control_parameters,
    contour_parameters=contour_parameters,
    electrode_calculators=
        [left_electrode_calculator, right_electrode_calculator],
    )

device_configuration.setCalculator(calculator)

device_configuration = nlread("Au-C6H4S2-Au.nc")[1]

calculator=calculator(
        electrode_voltages=(0.0*Volt, 0.0*Volt))

device_configuration.setCalculator(
      calculator(),
      initial_state=device_configuration)
device_configuration.update()
nlsave("Au-C6H4S2-Au-0.0.nc", device_configuration)

Here, the device_configuration.update() only take 2 steps to converge since it based on an already converged *nc file Au-C6H4S2-Au.nc that is generated by nlsave as

Code

device_configuration.update()
nlsave('Au-C6H4S2-Au.nc', device_configuration)

From this test, I guess the *nc file generated by nlsave  can be used as an initial guess and a restart checkpoint for the new calculation.   With best regards, Guangping

[Anders Blom]

That was a bit convoluted, but ... yes.

[zhangguangping]

Dear Anders. About the restart feature, I have been puzzled according to my test on this feature. http://www.quantumwise.com/publications/tutorials/item/502-restarting-stopped-calculations As the tutorial mentioned, I reivsed my script of input file to

Code

#----------------------------------------
# Device Calculator
#----------------------------------------
calculator = DeviceLCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=device_numerical_accuracy_parameters,
    iteration_control_parameters=device_iteration_control_parameters,
    contour_parameters=contour_parameters,
    electrode_calculators=
        [left_electrode_calculator, right_electrode_calculator],
    )

device_configuration.setCalculator(calculator)

device_configuration = nlread("checkpoint-0.4.nc")[0]

device_configuration.update(force_restart=True)

nlsave("Au-C6H4S2-Au-restart.nc", device_configuration)

transmission_spectrum = TransmissionSpectrum(
    configuration=device_configuration,
    energies=numpy.linspace(-2,2,100)*eV,
    kpoints=MonkhorstPackGrid(6,6),
    energy_zero_parameter=AverageFermiLevel,
    infinitesimal=1.36057e-05*eV,
    self_energy_calculator=RecursionSelfEnergy(),
    )
nlsave("Au-C6H4S2-Au-restart.nc", transmission_spectrum)
nlprint(transmission_spectrum)
transmission_spectrum.current()

The file checkpoint-0.4.nc is generated by CheckpointHandler in a previous run (and this job terminated normally, that is the SCF convergence has been done). Here, whether I use a checkpoint *nc file or a *nc file generated by nlsave, the calculation always does from the scratch. What wrong with my use of force_restart=True ? In constrast, for the following script,

Code

#----------------------------------------
# Device Calculator
#----------------------------------------
calculator = DeviceLCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=device_numerical_accuracy_parameters,
    iteration_control_parameters=device_iteration_control_parameters,
    contour_parameters=contour_parameters,
    electrode_calculators=
        [left_electrode_calculator, right_electrode_calculator],
    )

device_configuration.setCalculator(calculator)

old_calculation = nlread("Au-C6H4S2-Au-0.4.nc", DeviceConfiguration)[0]

device_configuration.setCalculator(
    calculator(electrode_voltages=(0.2*Volt, -0.2*Volt)),
    initial_state=old_calculation,
)
device_configuration.update()

nlsave("Au-C6H4S2-Au-initial2.nc", device_configuration)

transmission_spectrum = TransmissionSpectrum(
    configuration=device_configuration,
    energies=numpy.linspace(-2,2,400)*eV,
    kpoints=MonkhorstPackGrid(6,6),
    energy_zero_parameter=AverageFermiLevel,
    infinitesimal=1.36057e-05*eV,
    self_energy_calculator=RecursionSelfEnergy(),
    )
nlsave("Au-C6H4S2-Au-initial2.nc", transmission_spectrum)
nlprint(transmission_spectrum)
transmission_spectrum.current()

The sentance

Code

old_calculation = nlread("Au-C6H4S2-Au-0.4.nc", DeviceConfiguration)[0]

also can be written as

Code

old_calculation = nlread("Au-C6H4S2-Au-0.4.nc")[0]

and Au-C6H4S2-Au-0.4.nc file can be either the one generated by CheckpointHandler or nlsave, the continue calculation is OK. So, what is the story for the restart calculation by using force_restart=True? And what is the difference beween the *nc files generated by CheckpointHandler and nlsave for a continue scf calcualtion? With best regards, /Guangping

[Jess Wellendorff]

Dear Guangping, The option "force_restart=True" forces the self-consistent calculation to restart from scratch, even though the calculation may already be converged. Here's how I read the first script in your post below: 0: a device configuration was probably defined earlier in the script. 1: you define a device calculator. 2: you attach the calculator to the configuration. 3: you then re-define the device configuration by "nlreading" a saved nc file,  so the calculator defined in step 1 is no longer used (it was attached to the original configuration, which you just replaced). 4: you call "device_configuration.update(force_restart=True)", forcing the calculator to redo the SCF from scratch using whatever calculator that was saved in the NC file along with the configuration. 5: the configuration is saved, including the re-converged calculator, and you do the transmission spectrum. If an NC file contains a converged calculation, I suggest you just nlread that file and head straight for the post-SCF analysis:

Code

device_configuration = nlread("file.nc")[-1]
transmission_spectrum = TransmissionSpectrum(
    configuration=device_configuration,
    ....,
    )
nlsave("file.nc", transmission_spectrum)
nlprint(transmission_spectrum)
print transmission_spectrum.current()

[zhangguangping]

Dear Jess,

Thanks for your kind reply.

Please find the atttached the full *py input file. I just test the restart feature as I learnt from this tutorial http://www.quantumwise.com/publications/tutorials/item/502-restarting-stopped-calculations.
In this tutorial, it said that force_restart=True can be used to do the calculation from where it stoped not from the scratch.

As I have tested, if I use force_restart=True, I always to start the calculation from the scratch. And this is my puzzling.

So, there is something wrong in the tutorial, Right?

With best regards,

/Guangping