Tuesday 10 June 2014

5. Performance considerations when using nwchem

nwchem supports parallel computing among all the available CPUs in a multicore computer, depending on the compilation flags. The Ubuntu distribution is linked against openmpi library, but compilation flags also determine whether they benefit automatically from the available CPUs or whether you need to explicitly indicate how many CPUs to use. In order to check how your installation works, first check whether it depends on openmpi:
linux$ ldd /usr/bin/nwchem | grep libmpi
If you get a match, then you can use mpi parallelization. Next step is to check how many CPUs your computer has with lscpu command. and whether nwchem distributes the process among all available available ones. Run the command again and in another shell run top. If the number of CPUs is bigger than 1, but there is only one nwchem process, then you must explicitly run nwchem in multiple cores. You should run instead:
linux$ mpirun -n <number of cpus> nwchem l-proline
Not the entire optimization process is completely parallelizable, so parallelization does not meant that if you run the process in 4 cpus, it will take 4 times less time, but for sure will improve it. Improvement also depends on the kind of calculation you are performing. Some tasks are more paralellizable than others.

Once you achieved the optimized geometry, the output files are what you need for performing other calculations with it. But if you want as output a new xyz file for visualization or other purposes (even using this optimized geometry in a geometry block of another nwchem file), you can print it, by including

driver
  xyz l-proline_opt
end
 This command will generate a new xyz file for each optimization step performed. So, lets modify our nwchem input file:
basis
  * library 6-31g
end
driver
  xyz l-proline_opt
end
task scf optimize
(we don't need the geometry block here, as we will work with the generated vector files, but if you do include it, it will be ignored)

As we are working with an already converged solution, nwchem will perform some few more steps (and then discover that convergence was already achieved), so you will get some few l-proline_opt* files, consecutively numbered. The last one contains the optimized geometry. You can then do:
linux$ babel l-proline_opt-002.xyz l-proline.mol
and open the file with jmol for visualization.

You can also use the optimized geometry for a further task instead of the nwchem output vector files. You can generate the last exactly in the same way we generated them the first time, but with the optimized xyz file instead in the geometry block. As geometry is already optimized, nwchem will achieve convergence very fast, in a pair of steps. Even faster, you can use task scf energy instead of task scf optimize. This will generate all the vector files without performing any optimization.

A typical strategy for optimizing molecule geometry with a higher order (more accurate and more CPU demanding) basis set, is to pre optimize molecule geometry with a fast, less accurate one, so the optimization with the higher one will take really much less time. For example, after we got an optimized state with 6-31g, we can run a further one with cc-pvdz:
basis
  * library cc-pvdz
end
task scf optimize 

4. Energy optimization of molecule geometry with nwchem

In this lesson I will introduce a great opensource suite for classical and quantic chemical computation, called nwchem. If you are not familiarized with this application, I recommend to check the Getting Started section of the documentation before continue reading this post. Here i will mention only the most basic concepts.

nwchem reads all the parameters needed for the computation from an input file with extension nw. In order to feed nwchem with the molecule geometry, we can write the coordinates in the input file, in simple xyz format, or load it from another file. As in out first post, Visualizing a 3D structure of a molecule with jmol, we will start geometry specification from a smiles file. Lets suppose we want to minimize the configuration of L-proline (smiles: [O-]C(=O)C1CCC[NH2+]1), so we save the smiles string into a file, and convert to xyz using babel:
linux$ babel l-proline.smi l-proline.xyz --gen3d
And the following is the nwchem input file for optimizing configuration:
geometry
  load format xyz l-proline.xyz
end
basis
  * library 6-31g
end
task scf optimize
Inside the geometry directive we instruct nwchem to load molecular geometry from our xyz file.

Inside the basis directive we specify the basis set to use for modeling the electronic orbitals of each atom class in the molecule. There are many different basis sets. You have yo choose in terms of compromise between speed and precision. The set cc-pvdz is an acronym for the correlation-consistent Double-zeta, and it is the minimal of the correlation-consistent group of basis sets, so it is the faster and less accurate of them, but still good quality for most serious applications that involves atoms up to the 3rd period. However, in many systems still can demand lot of time, specially if the input xyz file is very far from optimized configuration. So for this exercise we will choose a quite fast (but less accurate) model: 6-31g. We apply this model to all the atoms, using the wildcard *.

With the task directive we invoke the Hartree-Fock (SCF) method for optimizing the molecular geometry by energy minimization.

Save the file as l-proline.nw and then run
linux$ nwchem l-proline
The computations are very CPU intensive. They imply quantum mechanical numeric methods with several levels of iterations, and increases geometrically with the number of atoms involved, so they will take a time, which depends mainly on the CPU speed on your computer and the geometry of the input file.

If you get an error saying that optimization failed because computation didn't reach convergence criteria, don't worry. Just run the command again, until you achieve the default set convergence threshold. nwchem stores iteration vectors state so when you run the command again, you will not start from beginning, but from the point the previous run finished. Even more, the saved vector files will allow you to perform further calculations (see following lessons).