QMMM: nwchem - options

[mbarbry]

It seems that to be able to use nwchem with QMMM, the option

set bq background

must be present. However, if the option is missing, no error messages are printed and the program get stuck into a while loop (it seems). At least nwchem calculations fail since the molecule geometry is not written in the input file.

A line checking if _write_charges is true would be necessary.

[JensWehner]

which branch?

[mbarbry]

I'm in master branch. The problem come from the line

if (_write_charges) {
     qmatoms = orbitals_guess->QMAtoms();
} else {
    QMMInterface qmmface;
    qmatoms = qmmface.Convert(segments);
}

in qmpackages/nwchem.cc

[JensWehner]

It is fixed in tutorial-fixes

[mbarbry]

ok, I will give it a try!

[mbarbry]

With the branch tutorial-fixes, there is still some error message missing. In my test, somehow, there is no background charge (probably caused by some issue with previous jobs), then nwchem fails with the error

bq_input: max_nbq                   0

However, this is still not visible from the xtp side which continues to run.

[JensWehner]

did it write the set bq background option

[mbarbry]

It writes some option related to bp, however set bq background is not present. system.nw.txt

[JensWehner]

okay let me check, give me 30 min

[JensWehner]

Did you have it in there before?

[mbarbry]

The options? In the master branch there were written only if I added set bq background to the options manually in dftpackage.xml, otherwise, nothing was written, not even the geometry of the molecule.

[JensWehner]

Yeah, I tried to fix it in tutorials_fix if you set set bq background manually there it should work just fine.

[mbarbry]

In my particular case, it does not work because the total number of charge is null. Thus nwchem failed to run since it expects some charges

0:bq_input: max_nbq:Received an Error in Communication

I think then, that an error should be return if the number of charges is null. In qmpackages/nwchem.cc we have the line

int numberofcharges=WriteBackgroundCharges(_crg_file);

it is then easy to add a check,

if (numberofcharges == 0.0) {
    die("no charge in the background, check your inputs, or whatever other message")
}

By the way, any clue why I don't have charges in my test? It is difficult to understand why from the code.

[JensWehner]

Good question. Are your cutoffs to small?

[JensWehner]

Ahh okay. I will add this exception to the code

[mbarbry]

No, I took quite a large cutoff and saw that quite a lot of molecules were participating to the background charge. I think it is because of my ab initio calculations. When looking to the mps file, I see no charge. segMethane_e.mps.txt

[mbarbry]

I think the partial charge reader for nwchem is outdated or I did not run my ab initio calculations correctly and I don't get the line expected by the parser ....

[JensWehner]

yes the mps file has no charge, so nothing is written out.

It should have an error. How did you make the mps file?

[mbarbry]

I used log2mps to generate the mps file. but I found out why It did not read the charges, it is because log2mps is expecting the string EPS in the nwchem log. If this string is absent then no charges are read. I have to modify thus my nwchem calculations in order that they perform electrostatic potential calculations.

[JensWehner]

Yes, log2mps cannot read charges without charges being calculated.

[JensWehner]

This is really valuable for us. You have at least discovered 3 edge cases so far.

[mbarbry]

Ja, I guessed that having a someone testing the code is nice for you to improve votca, for me it is a painful path :D

[JensWehner]

sorry. I do apologize.

[mbarbry]

No issues, I know what is it do develop scientific software

[baumeier]

Ja, I guessed that having a someone testing the code is nice for you to improve votca, for me it is a painful path :D

@mbarbry Just chiming in to say that it is really helpful. Especially with the ways the different qmpackages work, and nwchem not being our "work horse", there are a bunch of things we don't easily spot.

[mbarbry]

I'm glad to help! I found another bug. I'm getting a segmentation fault in gwbse/sigma.cc, routine Sigma::X_diag.

It is concerning the loop

for (unsigned i_occ = 0; i_occ <= _homo; i_occ++)

In my case, _homo = - 1 (probably because I'm not initializing my jobs correctly). Thus, the loop will never finish, and we finish to get a segfault when i_occ is too large.

You should probably had a test there,

if (_homo < 0) die("Error message")

[mbarbry]

or just using an int instead of an unsigned for i_occ would avoid entering into the loop. But my knowledge of C++ is rather limited, you probably know better how to handle such situation.

[mbarbry]

The issue is that in this particular calculation, the number of electrons detected by xtp is 0. This is because of the orbitals file converted with the mov2asc executable. The number of electrons is obtained by summing over when a number in this file is equal to 2. qmpackages/nwchem.cc

_input_file >> _occ[ _imo ];
if (_occ[ _imo ] == 2.0) {
    _number_of_electrons++;
}

Somehow in my file, I have 1 instead of 2, thus the number of electrons is null. Any clue why?

[mbarbry]

Do you know why the number of electron is read from the molecular orbitals file in nwchem? It is directly accessible from the main output.

No. of electrons :    10
Alpha electrons :     5
Beta electrons :     5

system.log

[JensWehner]

Yes, one flaw of votca is that is was originally designed for closed shell molecules. E.g. alpha electrons are only counted. It is flawed

[mbarbry]

Ok, I just opened a pull request. I think there is a bug in the nwchem parser.

[baumeier]

What I find odd is that this part has always worked for closed-shell systems since 5 years ago for the calculation of transfer integrals. Which version of nwchem does this happen with? Maybe they changed the structure of the file, putting occ=1 and having alpha and beta electrons separately.

Here one of my outputs (ignore the nitrogen blabla comment, it's been a different system) from NWchem 6.6:

`# This is an NWChem movecs file translated by mov2asc a4d7e69728df14518955b70354ce20bb cf77caddb7da017694d26f313dcfa196

dft Tue Apr 24 19:32:33 2018 dft 42 Nitrogen cc-pvdz SCF geometry optimization 8 ao basis 1 535 535 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.200000000000000E+01 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00`

[mbarbry]

The thing is that to do calculations with charges, we need to do open shell calculations? Or maybe I did not understand correctly how to do calculations with charges with nwchem.

[baumeier]

Do you want to do GW on a charged system? That is currently not implemented, but on our list.

[mbarbry]

No, not for the moment at least. This may explain why it is not working. I was getting an error running GW with nwchem, I was thinking it was because it was running DFT with charged systems. I see a folder created by xtp xjob_2_1:segmethane:e which run DFT with a charge. Then, to perform this calculation, nwchem uses openshell, and then the nwchem parser does not manage to gets the number of electrons.

system.nw.txt

[baumeier]

Yeah, if you run GW on open-shell systems, this is sort-of expected behavior. In the sense that it should fail, ideally with a check before.

Now, why it is doing this inside xjob_2_1:segmethane:e at all with GW, is something I don't really understand. We can run plain DFT in QMMM with nwchem also for charged molecules but it relies on nwchem's internal espfit procedure. Then the parser should only read energy and esp charges, or at least only use those.

[baumeier]

OK, you must have jobs in the qmmm.jobs file with segmethane:e and segmethane:h or something, while at the same time a section in qmmm.xml options files. That's at least how I can rationalize that this combination can occur.

[JensWehner]

The problem is probably a mismatch between the tags in qmmm and gwbse

[mbarbry]

A solution to avoid issues in reading the number of electron with nwchem would be to get them directly from the main output instead of the molecular orbitals. They are written in plain text so no rooms for mistakes.

No. of atoms     :     5
No. of electrons :    11
Alpha electrons :     6
Beta electrons :     5
Charge           :    -1
Spin multiplicity:     2

[baumeier]

That would seem a sensible option to me. Plus for the moment that if the parser tries to read MOs with occ=1, throw an error - until we have a data structure to read alpha and beta orbitals into separately.

[JensWehner]

I close this issue, as it basically pertains to the larger issue of spin. THe background charge problem was fixed