561. b3pw91 in nwchem and g09

UPDATE: there was an error in the earlier version where I gave the wrong energy for the b3pw91 functional in nwchem. In the old version the energy I provided was very close to that of acm in nwchem rather than b3pw91 in g09.

Note that for a large molecule with a medium sized basis set (101 atoms, ca 1100 functions,  ca 2200 primitives) the energy difference between b3pw91 in g09 and b3pw91 in nwchem as defined below is 0.0124 Hartree, which is pretty big (7.8 kcal/mol), although in absolute terms it's quite small (nwchem: -6187.741840960054 Hartree. g09: -6187.75427966 Hartree).

The difference is a lot smaller for the small molecule in the example below.

Original post:
According to http://www.nwchem-sw.org/index.php/Special:AWCforum/st/id721/Are_these_definitions_correct_fo....html b3pw91 (as defined in Gaussian 09) and acm (as defined in nwchem) are identical.

Looking at the energies I've been getting, that's not true when it comes to G09 and NWCHEM.

---

That acm and b3pw91 are the same should be reasonable -- b3 indicates that it's Becke's 3-parameter hybrid exchange correlation functional model, which is also known as the Adiabatic Connection Method (ACM).

For historical reasons, g98 implemented the ACM as B3LYP, by using LYP instead of PW91, and using VWN_1_RPA and a few other tricks -- see section 2 in http://verahill.blogspot.com.au/2013/06/446-b3lyp-and-wah-confusion.html

Then it would stand to reason that B3PW91 would be the 'canonical' version of Becke's 3-parameter functional.

Looking at http://www.nwchem-sw.org/index.php/Release62:Density_Functional_Theory_for_Molecules acm is defined as

xc HFexch 0.2 slater 0.8 becke88 nonlocal 0.72 vwn_5 1 Perdew91 0.81

(there are several versions of VWN -- I know it's vwn_5 from the output)

---

Either way, using acm in a single energy calculation (no optimisation) in nwchem on a water molecule with 6-31+G* (acm/6-31+G*) gives

-76.358375905073 Hartree

G09 using B3PW91/6-31+G* (manually defined basis set so we're using the same form in both nwchem and g09) gives

 -76.3557851653 Hartree

nwchem using

xc HFexch 0.2 slater 0.8 becke88 nonlocal 0.72 vwn_5 1 Perdew91 0.81

gives

-76.358375905072 Hartree

and nwchem using

XC HFexch 0.20 slater 0.80 becke88 nonlocal 0.72 perdew91 0.81 pw91lda 1.00

obtained from http://myweb.liu.edu/~nmatsuna/gamess/refs/howto.dft.html,gives

-76.355784373093 Hartree

This last definition is thus equivalent to b3pw91 in g09.

The gaussian manual is less than helpful. In fact it is quite misleading:

"These functionals have the form devised by Becke in 1993 [Becke93a]:
A*E_X^Slater+(1-A)*E_X^HF+B*ΔE_X^Becke+E_C^VWN+C*ΔE_C^non-local
[..] B3LYP uses the non-local correlation provided by the LYP expression, and VWN functional III for local correlation (not functional V). [..]B3P86 specifies the same functional with the non-local correlation provided by Perdew 86, and B3PW91 specifies this functional with the non-local correlation provided by Perdew/Wang 91.

---

Addendum:
While I think B3PW91 should be the same as ACM in nwchem (note that nwchem does not have b3pw91 as a keyword), I decided to have a look at how different packages define b3pw91.

nwchem -- doesn't exist. Manual.

g09 (this post) -- xc HFexch 0.20 slater 0.80 becke88 nonlocal 0.72 perdew91 0.81 pw91lda 1.00

gamess US (here) -- xc HFexch 0.20 slater 0.80 becke88 nonlocal 0.72 perdew91 0.81 pw91lda 1.00

PQS (page 52, manual) Paraphrased:
"B3PW91 -- hybrid 3-parameter HF-DFT functional comprising combination of Slater local exchange, Becke nonlocal exchange, VWN 5 local correlation and PW91 nonlocal correlation together with a portion (20%) of the exact Hartree-Fock exchange (original 3-parameter hybrid recommended by Becke)". That to me sounds like ACM.

Turbomol -- not available. Manual.

Orca -- "B3PW The three-parameter hybrid version of PW91". Not informative.

molpro -- doesn't exist. manual

Dalton -- (page 285, manual).

"B3PW91 3-parameter Becke-PW91 functional, with PW91 correlation functional. Note that PW91c includes PW92c local correlation, thus only excess PW92c local correlation is required (coe cient of 0.19).
Combine HF=0.2 Slater=0.8 Becke=0.72 PW91c=0.81 PW92c=0.19"

So the local correlation is 1*PW92c= 0.81 PW91c + 0.19 PW92c. This is, I presume, is quite different from VWN.

Q-Chem -- "B3PW91 (B3 Exchange + PW91 correlation)". Not explicit enough for me.

560. Globus and XSEDE/Stampede on linux

I'm currently applying for extra computing resources via XSEDE through my US collaborators.

While joining XSEDE was as easy as signing up to any other website, and the system with a Campus Champion works really well, transferring files to and from your desktop and the compute resource (here stampede.tacc) is a bit trickier.

As an XSEDE user you can't transfer files directly to the compute resource via sftp.

I first tried the XSEDEFileManager.jar XUP application, and while it connected fine to stampede, it wouldn't allow the transfer of any files. Maybe it has to do with me using openjdk, maybe not.

I then tried the browser version of XUP, but it wouldn't even allow me to connect to a compute resource.

I then tried Globus Online, which is a service for transferring large data files which is managed by the University of Chicago.

1.The first step is to sign up for Globus Online. You'll get a confirmation email, click a link and you're up and running.

2. Follow this: https://support.globus.org/entries/23881557-Globus-Connect-Personal-for-Linux.
Install the desktop client so that you can set up your desktop as an endpoint. That involves several steps:

A. make sure you have tcl/tk installed. On Debian

sudo apt-get install tk tcl tcllib

B. make sure that you have ssh keys

ls ~/.ssh/id*

There should be a file/files called id_rsa and/or id_dsa as well as the same files with .pub endings.

If you don't, then generate ssh keys using
ssh-keygen

C. Log on to globus online and go to manage endpoints, click on "add Globus Connect Personal".

Enter a name for the system, e.g. Beryllium

A setup key will be generated -- you'll need this to set up the desktop client.

Click on the Linux link to download the linux client (globusconnect-latest.tgz)

D. In e.g. your ~/Downloads, do

tar xvf globusconnect-latest.tgz
cd globusconnectpersonal-2.0.2/
sudo mv globusconnectpersonal-2.0.2 /usr/local/bin
/usr/local/bin/globusconnectpersonal-2.0.2/globusconnectpersonal

Enter the code. You can then reload the end point manager window and the status should be changed to ready.

3. Next you'll go to 'Transfer Files' on the Globus Online website. Pick the right endpoints and transfer away!

559. Briefly: Nudge Elastic Band in NWChem

I'm currently exploring a method called Nudge Elastic Band to find the minimum energy pathway in a reaction involving a large metal cluster. While the NWChem documentation isn't bad, it could be clearer for happy amateurs like myself, which is the impetus for this post.

As usual, this is how I understand it -- which may be wrong.

The NEB method is described here.

So...

There are a number of ways of modelling reaction pathways computationally.

Brute force PES scan (g09 or nwchem)
The fastest, cheapest, least accurate one would be to make a simple Potential Energy Surface (PES) scan by keeping all coordinates frozen with the exception of the ones directly involved in a reaction
e.g. for Me-OH -> Me + OH you'd simply vary the C-O distance.

If you don't let the structures relax for each change in bond length you will overestimate the activation energy. If you do let the structures relax you run the risk of every single snapshot relaxing back to either pure methanol or pure Me and OH, rather than giving you a series of transitional energies.

Transition state search e.g. QST2 (g09)
A somewhat more sophisticated method is to try to identify a proper transition state (TS). You can do this by making a lucky guess and try to do a transition state optimisation. You can improve your chances by using e.g. QST2 (or QST3, which also uses a specific TS guess), which takes the product and the starting material and interpolates the coordinates to generate a transition state guess. This has worked pretty well for me for cases where a simple transition state is expected (i.e. a hydride transfer). You can then model the reaction path by doing an IRC calc.

Minimum energy pathway methods (chain-of-states) e.g. Nudge Elastic Band (nwchem)

A more sophisticated method is to use a chain-of-state method such as Nudge Elastic Band (NEB). There's a nice page about it here: http://theory.cm.utexas.edu/henkelman/research/saddle/

NEB generates a series of structures based on the starting point and the product. These are initially generated a simple interpolations between the coordinates of the starting point and the end point -- this is very similar to a brute for PES scan and gives an unnaturally high transition energy. However, the individual structures are then optimised in passes.

10 passes. Structural convergence is not achieved.

Not clear? Well, here's a rough algorithm:

1. Generate a series of linearly interpolated structures (beads) based on the Starting Point (SP) and End Point (EP) structures. In the simplest case, for N structures, each structure n is SP+(n/N)*(EP-SP).

Now, typically you're more interested in structures in the middle than at the ends (i.e. to get better resolution for the transition) and that is manipulated using the spring constant, kbeads. See figure 3 in http://theory.cm.utexas.edu/henkelman/pubs/jonsson98_385.pdf for an example of varying the spring constant. 1.0 seems to be reasonably safe.

The number of structures is set by nbeads. Make sure to use a reasonable number so that you get enough resolution.

2.  Each structure thus generated is submitted to a single geometry optimisation step i.e. it's not optimised until convergence, but only takes a single step along the energy gradient. The speed at which it is decending along the energy gradient is set using the stepsize.

The largest possible stepsize will give the fastest descent (which is good), but too large a stepsize may cause issues due to the atoms being moved too far (which is bad) and in turn cause SCF convergence failures due to the structures becoming unreasonable. So try 1.0 and hope for the best.

Rule of thumb: if you decrease the stepsize by a factor of 10 from 1.0 to 0.1, you should increase maxiter by a factor of 10.

3. The energies of the new structures are calculated and a reaction profile is generated.

4. Repeat step 2 and 3 until structural convergence or maxiter steps.



5. If you're lucky your calculation will end due to structural convergence. If not, and if you want to restart your calculation you can read in the last set of structures using xyz_path.

But make sure that you do print the intermediate structures in the first place by using print_shift 1.