Re: About changing charge/multiplicity during QM-MM

From: Gerard Rowe (GerardR_at_usca.edu)
Date: Tue Apr 02 2019 - 11:14:57 CDT

Francesco,

I'm not sure I'd call 0.5 spin population negligible. I see spin population values of about that magnitude for Cu(II) centers in metal complexes. What I'd take a close look at is where that unpaired spin resides. It's very common for spin density to bleed out from a metal onto its immediate ligand atoms. If the unpaired spin density is a few atoms removed from the metal, it's possible that you're seeing a localized ligand radial. Depending on what you're trying to get out of these simulations, it could be worthwhile to search for a broken-symmetry solution to the wavefunction.

-Gerard

________________________________
From: Francesco Pietra <chiendarret_at_gmail.com>
Sent: Tuesday, April 2, 2019 10:36 AM
To: Gerard Rowe; NAMD
Subject: Fwd: namd-l: About changing charge/multiplicity during QM-MM

As far as it matters, I have to correct my mail: actually the redox partner carries negligible spin (0.5, against 3.9 for both iron ions), while, unbound, it was a triplet. Previously I thought that 0.5 might be significant, while probably not. At any event, all other atoms carry out <0.00 spin, except one O-atom of iron-coordinated GLU, which carries 0.23 spin.
fp

---------- Forwarded message ---------
From: Francesco Pietra <chiendarret_at_gmail.com<mailto:chiendarret_at_gmail.com>>
Date: Sun, Mar 31, 2019 at 5:55 PM
Subject: Re: namd-l: About changing charge/multiplicity during QM-MM
To: Gerard Rowe <GerardR_at_usca.edu<mailto:GerardR_at_usca.edu>>

Hi Gerard
It was a problem with the cluster. By setting

--nodes=1
--ntasks=1
--cpus-per-task=1

the singlepoint had no problem.

I can read the Mulliken (or Loewdin) spin populations, so that I know where the spin mostly resides (also the redox partner carries some spin), but I am unable to get a plot of spin densities. With orca_plot I could generate a cube file that is read by AVOGADRO but the output is limited to molecular orbitals; apparently it does understand the alpha-beta difference.. With VMD things went on even less satisfactorily. Any hint?
thanks
francesco

On Fri, Mar 29, 2019 at 5:09 PM Gerard Rowe <GerardR_at_usca.edu<mailto:GerardR_at_usca.edu>> wrote:
If it still looks to a GBW for an initial guess, then you should be able to circumvent that by specifying an alternative initial guess method. Try adding PMODEL to your standard input line.

________________________________
From: Francesco Pietra <chiendarret_at_gmail.com<mailto:chiendarret_at_gmail.com>>
Sent: Thursday, March 28, 2019 12:49 PM
To: Gerard Rowe
Cc: namd-l_at_ks.uiuc.edu<mailto:namd-l_at_ks.uiuc.edu>
Subject: Re: namd-l: About changing charge/multiplicity during QM-MM

Gerard
Do you know how to set the keyword "noautostart" to orca. The keyword is mentioned in the instructions about restarting, however I was unable to get it accepted. For example "%scf noautostart end" failed. The point is that orca invents a gbw file even if none is present, and crashes because it does not find what was expected.
thanks
francesco

On Tue, Mar 26, 2019 at 10:24 PM Gerard Rowe <GerardR_at_usca.edu<mailto:GerardR_at_usca.edu>> wrote:
Unless your QM region is prohibitively large, it might be easier to just perform a single point calculation from scratch with the coordinates in the input file (i.e., remove the GBW file). That should avoid any issues.

________________________________
From: Francesco Pietra <chiendarret_at_gmail.com<mailto:chiendarret_at_gmail.com>>
Sent: Tuesday, March 26, 2019 4:55 PM
To: Gerard Rowe
Subject: Re: namd-l: About changing charge/multiplicity during QM-MM

Gerard
With full path to orca (4.1.1, the same path user for namd qmmm) same problem about the executable. In addition, a (for me) confusing allusion to basis sets.
francesco

WARNING!!!!!!!
Your GBWFile is either corrupt or from a different ORCA version!
Please be VERY careful with your calculation results!!!

Data_size, sizeof( TOrcaInfo ): 497844, 491312
     ... Fine, the calculation information was read
     ... Fine, the file contains a basis set
[file orca_tools/Tool-GTO-Bases/basissets.cpp, line 703]: TBasisList: Unable to read number of basis sets!

On Tue, Mar 26, 2019 at 6:09 PM Gerard Rowe <GerardR_at_usca.edu<mailto:GerardR_at_usca.edu>> wrote:
Hmm, it kind of looks like you either have a different number of atoms in your input than the program expects, or there are two different versions of orca installed on your system. You could try specifying the full path of Orca when running the program. That's the recommended way to do it for oMPI purposes, anyways.

________________________________
From: Francesco Pietra <chiendarret_at_gmail.com<mailto:chiendarret_at_gmail.com>>
Sent: Tuesday, March 26, 2019 12:35 PM
To: Gerard Rowe
Cc: NAMD
Subject: Re: namd-l: About changing charge/multiplicity during QM-MM

Hi Gerard:
I copied the files that you indicated to a new folder. Then I modified qmmm_0.input by deleting "enGrad" from
the first line of the .input file "! UKS BP86 RI SV def2/J enGrad KDIIS SOSCF"
saving it as qmmm_0.noGrad.input, then launched a job

srun orca qmmm_0.noGrad.input > qmmm_0.noGrad.log

It crashed with message from .log file:

Checking for AutoStart:
The File: qmmm_0.noGrad.input.gbw exists (it was renamed by the system)
Trying to determine its content:

Checking for AutoStart:
The File: qmmm_0.noGrad.input.gbw exists
Trying to determine its content:
     ... Fine, the file contains calculation information
     ... Fine, the file contains calculation information
     ... Fine, the file contains calculation information
     ... Fine, the file contains calculation information

WARNING!!!!!!!
Your GBWFile is either corrupt or from a different ORCA version!
Please be VERY careful with your calculation results!!!

Data_size, sizeof( TOrcaInfo ): 744532, 491312
     ... Fine, the calculation information was read
     ... Fine, the file contains a basis set

WARNING!!!!!!!
Your GBWFile is either corrupt or from a different ORCA version!
Please be VERY careful with your calculation results!!!

Data_size, sizeof( TOrcaInfo ): 744532, 491312
     ... Fine, the calculation information was read
     ... Fine, the file contains a basis set
[file orca_tools/Tool-GTO-Bases/basissets.cpp, line 706]: Wrong number of basis-sets stored!

[file orca_tools/Tool-GTO-Bases/basissets.cpp, line 706]: Wrong number of basis-sets stored!

[file orca_tools/Tool-GTO-Bases/basissets.cpp, line 706]: Wrong number of basis-sets stored!

[file orca_tools/Tool-GTO-Bases/basissets.cpp, line 706]: Wrong number of basis-sets stored!
____
The .err file said I/O operation failed.
_________

I don't understand that "Wrong number of basis set". The namd-qmmm simulation had been interrupted by the walltime.
Could that be a reason for orca corrupted files?

thanks
francesco

On Mon, Mar 25, 2019 at 7:34 PM Gerard Rowe <GerardR_at_usca.edu<mailto:GerardR_at_usca.edu>> wrote:
Francesco,

It sounds like you may have a fortunate outcome. I don't know the nature of your redox partner, but it's not an unreasonable assumption that the spin state will be the same between reactant and product if your partner becomes a ligand. If the redox partner started as a radical and became a closed-shell molecule, your final wavefunction is probably going to behave, and the trajectory could be well-suited for analysis through molecular orbital plots or spin-density plots. I think the most useful information will come out of spin density for a redox reaction. Things get more complicated (but not impossible) if the redox partner has unpaired electrons in it.

As for carrying out stand-alone Orca calculations from the temp folder files of QM/MM, you should copy the inp, gbw, and point charge files into a new folder. Unless you have a lot of experience with QM programs, I'd stick to single point calculations here because you would have to play a tricky game of atomic restraints for QM geometry minimization.

Your choice of level of theory is pretty good for QM/MM, but if I were trying to get a more reliable wavefunction, I'd probably change the basis set to def2-svp, and possibly specify def2-tzvp for the metal ion. If you don't have restrictions on how much CPU time you use for QM/MM, you may also play around with other local DFT functionals like PBE or M06L to see if they lead to different equilibrium structures. M06L, in particular, is parameterized to handle non-covalent interactions.

-Gerard

________________________________
From: Francesco Pietra <chiendarret_at_gmail.com<mailto:chiendarret_at_gmail.com>>
Sent: Monday, March 25, 2019 12:09 PM
To: NAMD; Gerard Rowe
Subject: Re: namd-l: About changing charge/multiplicity during QM-MM

Gerard:
The redox partner is included but, regrettably, I have not yet carried out an analysis of what orca could tell. I simply observed that now the iron ions and the redox partner have reached a bonding distance, which remains stable on long QM-MM simulations, and the type of complex that is formed is experimentally known to have FeIII, never FeII. So, I came here much too early with my question, but, to go on, I'll take into account what you said. At any event the type of ligands of the Fe ions (in a protein) are for high spin, which over-complicates the matter. In addition, the simulation, because of the type and size of the system, was started at SlowConv but was continued at SOSCF in order to get convergence . All over it was at what might well be a too low level of DFT ("! UKS BP86 RI SV def2/J enGrad KDIIS SOSCF").
Thanks
francesco

On Mon, Mar 25, 2019 at 3:53 PM Gerard Rowe <GerardR_at_usca.edu<mailto:GerardR_at_usca.edu>> wrote:
Francesco,

If your metal ion is definitely undergoing redox change, that's going to be very difficult to treat with QM/MM unless you have also included the redox partner in the QM region. Even then, it's going to be problematic because the spin state of your resulting system is unlikely to be reasonable. This is especially true if your final system is supposed to have unpaired electrons in different regions that aren't in communication anymore. What kind of system are you studying?

I can imagine two scenarios where this will work out:

  1. Your Fe(II) state is low-spin d6 (singlet) and the redox partner is a doublet. In the product state, The Fe(III) is a low-spin doublet and the redox partner is now a singlet.
  2. Your Fe(II) state is high-spin d6 (sextet) and the redox partner is a singlet. In the product state, The Fe(III) is a high-spin quintet and the redox partner is now a doublet.

If scenario 1 applies, you're in luck, and the spin system will be reasonable. If it's scenario 2, you're going to run into serious issues unless the redox partner happens to be bound to the metal. Long-range separation of unpaired spin in a DFT calculation often leads to problems and non-physical solutions. You can check the spin density of your resulting wavefunction in a similar way that you can visualize molecular orbitals.

-Gerard

________________________________
From: owner-namd-l_at_ks.uiuc.edu<mailto:owner-namd-l_at_ks.uiuc.edu> <owner-namd-l_at_ks.uiuc.edu<mailto:owner-namd-l_at_ks.uiuc.edu>> on behalf of Francesco Pietra <chiendarret_at_gmail.com<mailto:chiendarret_at_gmail.com>>
Sent: Monday, March 25, 2019 3:16 AM
To: NAMD
Subject: namd-l: About changing charge/multiplicity during QM-MM

I dare posing a general question on QM-MM, surely a naive question to those experienced with QM-MM (I am doing such simulations for the first time).

I am observing a change in my metalloprotein system where apparently the iron ions become bounded to other ligands, changing from FeII to FeIII. The overall charge/multiplicity is not maintained, as far as I can understand. Nonetheless, ORCA is repeating its SCF cycles, converging after a number of cycles. I have already restarted the simulation six times for 24hr on 36 cores with the same settings, in particular as to the charge and multiplicity, and the structure of the complex is no more changing, as far as it can be judged from geometry. I can't figure out how the system readjusts itself so that the QM calculation can go on.
Thanks for any comment
francesco pietra

This archive was generated by hypermail 2.1.6 : Thu Dec 31 2020 - 23:17:10 CST