The ej_calc code, is divided into 4 modules:

1. main.F90: controls the flow of the code.

2. init.F90: declares all the global variables, reads the input files and initializes the global variable using the information obtained from the input files.

3. modul_lib.F90: builds up the Hamiltonian based on the definition obtained from the input file. It further forms different combinations of the equations and solves them using LAPACK routines. The average of all the valid sets of j-values and the standard deviations for each J-value is also determined using this module.

4. enrg_frm_jval.F90: calculates the energy of the all spin states using the Heisenberg Dirac van Vleck Hamiltonian (Ĥ = -2 ∑_i>j J_ijs_is_j). [Heisenberg], [Dirac], [Vanvleck]

A detailed breakdown of each module is provided below.

This module defines all the global variables.

Variable	Function
num_mag_cent	stores the number of metal centres in each configuration.
no_of_j_val	stores the total number of j-values defined for the system.
tot_interac	stores the number of possible interactions between spins i.e. the total number of S1S2 terms.
num_spin_dens_set	stores the number of sets of spin density available which is the total number of lines in the spin density file.
try	this helps in keeping the program running just in case there is an equation set is encountered that cannot be solved.
poss_comb	stores the possible combination of equations i.e. the (n-1)Cj term.
main_cntr	keep track of total non-singular equations.
start, finish	These in combination determine the time taken during the program execution.

Array	Function
hamil1, hamil2	These in combination keep track of which centres are interacting .
jpos	keeps track of how many interactions are being included under a single j-value.
energy	stores the energy for each electronic configuration.
tot_avg	stores the average of all the J-value sets which are valid and which meet the standard deviation criteria.
stddev	stores the standard deviation for the different J-values.
std_per	stores the percentage standard deviation for the different J-values.
spin	stores the spin operator value associated with each metal centre as given in the spin operator file.
hamil	stores the 2S1S2 type weight terms for each interaction.
jval	stores the total (2S1S2 type) coefficient associated with each j-value.
jval_trkr	stores all the computed sets of j-values.
singul	(for ej_calc_spin only) stores the location of singular sets of equations identified by ej_calc_form. It essentially stores the data in the file singul.

This module contains all the subroutines required to calculate all the possible sets of j-values and determine the valid ones based on standard deviation. The following subroutines form part of this module:

1. jvals: This subroutine will determine the S1S2 terms and which of these terms comes under a given j-value.

2. solv: This subroutine forms all possible sets of equations and solves them. The local variables used in this subroutine are as follows:

Array	Function
ctr	allows the separation and hence the evaluation of all possible combinations.
loc_energy	stores the energy locally since that will be changed in each set of equations considered.
enrg_tb_slvd1	local copy of the global array ‘energy’.
enrg_tb_slvd	stores the energy of the set of equations that will be solved in one instance.
enrg_tb_slvd2	(for ej_calc_form only) copy of the array ‘enrg_tb_slvd’.
loc_jval	store the coefficients of j-values locally since that will be changed in each set of equations considered.
jval_tb_solvd1	local copy of the global array ‘jval’.
jval_tb_solvd	stores the coefficients of j-values of the set of equations that will be solved in one instance.
jval_tb_solvd2	(for ej_calc_form only) copy of the array ‘jval_tb_solvd’.

In this subroutine, first of all, an electronic state is chosen as the reference state (every state is chosen as a reference state once) and a new set of equations is formed by subtracting the reference state from all the other states. These states are then grouped into all possible sets of r (the number of j-values one is looking for) and each such set is solved as simultaneous linear equations using the LAPACK library routine dgesv. The LAPACK routine dgesv is implemented in a different subroutine called simul_solv which also performs a few basic checks for singularity. If the solution is found to be singular, the value of all j-values for the given set is set to zero. If, however, the set is not found to be singular, a more thorough screening of the set is performed using the LAPACK library routine dgesvx which is implemented in the subroutine simul_solv1. Additionally, if any of the j-value is found to be larger than 50 cm^-1 for a given set of equations, details about the set and the corresponding J-values obtained are printed in the file bigg. In case of ej_calc_spin however, only those sets that are determined as non-singular by the array singul (generated from the 'singul' file created by executing ej_calc_formal) are solved.

3. checkpos: This subroutine is a helper subroutine for the subroutine ‘solv’ to cycle through all the possible combinations of r (the number of j-values one is looking for) spin states from n (or more accurately n-1) spin states.

4. simul_solv: This subroutine is also a helper subroutine for the subroutine ‘solv’. It checks if there are any identical pair of equations which will automatically make the whole set singular. If such pairs are not found, then a solution is calculated using the LAPACK library routine dgesv.

5. simul_solv1: (for ej_calc_form only) This subroutine is also a helper subroutine for the subroutine ‘solv’. It calculates the solution for the same set of equations that were considered solvable by simul_solv using the LAPACK library routine dgesvx which looks further for sets of equations close to singularity. The reason why this is used only for the sets which are found to be solvable by dgesv because it is more computationally expensive (The speed gain in keeping simul_solv and simul_solv1 has not been determined and considering that this code can still solve ~1,000,000 sets of equations within minutes, it was not considered worth doing).

6. soln_chk: (for ej_calc_form only) This subroutine makes sure that the code does not get stuck on any set of equations.

7. avrg_calcs: This subroutine calculates the average of all the valid set of solutions.

8. std_dev: This subroutine calculates the standard deviation for all the valid set of solutions.

9. std_dev1: This subroutine checks if all the valid j-value sets are within 3 standard deviations or not and if they are not, they are removed.

10. avrg_calcs1: This subroutine will calculate the new average for j-values since some have been eliminated by the subroutine ‘std_dev1’.

11. final_print: This subroutine prints the final set of J-values.

12. backtrack1: This subroutine will calculate the energies using the computed J-values and compares them with the energies that are provided as input.

This module calculates the energy of the spin states using the coupling constants calculated by modul_lib.F90. The main local variables used in this module are as follows:

Array	Function
unp_elec	stores the spin operator value for each metal centre.
all_config	stores all the possible spin configurations of spin states.
new_nrg	stores the energy of all possible spin states.
Variable	Function
poss_comb1	stores the total possible spin states.

The following subroutines form part of this module:

1. avg_spin: This subroutine calculates the average spin operator (in this case, the formal spin) value for each metal centre.

2. all_poss_spins: This subroutine along with the subroutine ‘update_spin’ forms all possible spin configurations, stores them in the array ‘all_config’, determines the energy of each configuration in accordance with the specified hamiltonian and stores the energy in the array ‘new_nrg’.

3. sort: This subroutine sorts all the spin states on the basis of increasing energy using bubble sort.

4. swap: This is a helper subroutine for the subroutine ‘sort’.

5. prnt: This subroutine prints the sorted energies in the nrg_frm_jval file.

6. update_spin: This is a helper subroutine for the subroutine ‘all_poss_spins’ for cycling through all possible spin configurations.

7. std_dev2d: (for ej_calc_spin only) calculates the standard deviation for the spin operator values that are used to calculate the energy of all spin states.