LOOS for a PCA - LAMMPS output

[gustavmf1]

Hi,
I’m trying to use LOOS for a PCA (with --svd) after an MD simulation in LAMMPS.
The MD simulation output should be in a DCD or PDB format, correct?
https://github.com/GrossfieldLab/loos/wiki/Using-LOOS-for-Principal-Component-Analysis
In a python script I should run:

import loos
--svd --noalign TrajectoryFileName.CorrectFormat

Is this the idea?

Thanks!
Best,
Gustavo

[agrossfield]

Hi,

To do PCA using loos, you don't need to write a python script -- there's an existing program (svd) that will do it for you. Re-reading the wiki page, I realize that's not as clear as it should have been.

You'll need to build and install loos, and as you said working from lammps you'll need a pdb file as your system file and a dcd as your trajectory.

The simplest command line would be svd foo.pdb bar.dcd, which would perform PCA on the CA atoms in foo.pdb, using the trajectory bar.dcd.

To see some example command lines and get a detailed explanation of the output, run svd --fullhelp

You'll want to be very careful about specifying --noalign: the usual use of PCA assumes that translation and rotation have been removed (producing 6 modes with 0 eigenvalues). svd ensures that by performing an iterative alignment. The --noalign option is there in case you know you're going to do a bunch of svd calculations, in which case you might want to pre-do the alignment instead of doing it on the fly. You shouldn't run svd on an unaligned trajectory unless you know exactly what you're doing -- the odds are it's not what you want.

Hope this helps,

Alan

[gustavmf1]

The system file is in the .xyz format:

9370
Atoms. Timestep: 0
1 18.8009 45.6819 209.906
2 19.472 44.3034 208.758
1 19.6258 44.2128 206.756
2 21.3715 44.6588 205.734
1 20.1638 44.7336 201.209
2 21.4423 44.9688 199.756
1 21.674 45.4212 196.621
2 21.6631 45.449 194.667

How can I change its atoms to atoms named CA? Do you have any suggestion or a link that shows how to proceed?
Also, is it necessary to specify the aligning in a command?

Thank very much for all the help! It hasn't been easy to find information.
Best,
Gustavo

[agrossfield]

Like I said in my previous post, when LOOS reads XYZ it's looking for Tinker format, not the format used by a qm package. What you're using is the XYZ format from a QM package, which loos can't read.

So, your first step is to use a file format that loos actually supports -- I gave a link in my previous post to what loos supports. For example, you can probably get lammps to convert to a PDB for you, and presumably that PDB will have atom names in it.

Second, you'll need to work out a selection for the atoms you're interested in. As a rule, you don't want PCA on all of the atoms, just a subset (eg PCA on solvent atoms is generally going to be meaningless). If you don't have atom names to work with, then you'll have to figure out which atoms are interesting, identify them by number, and make a selection string out of the atom indexes. If you do have atom names, then for a protein the default is 'name == "CA"', although you could plausibly want the selection 'backbone', which works for proteins and nucleic acids.

For more on the loos selection syntax, see http://grossfieldlab.github.io/loos/selections.html

[gustavmf1]

Hi Alan,

We've manage to create a .pdb from a .xyz file
My .pdb file looks like:
CRYST1 0.000 0.000 0.000 90.00 90.00 90.00 P 1 1
ATOM 1 O X 1 23.024 49.373 211.635 0.00 0.00 O
ATOM 2 ZN X 1 23.024 49.373 209.622 0.00 0.00 ZN
ATOM 3 O X 1 23.024 49.373 206.328 0.00 0.00 O
ATOM 4 ZN X 1 23.024 49.373 204.316 0.00 0.00 ZN
ATOM 5 O X 1 23.024 49.373 201.021 0.00 0.00 O
ATOM 6 ZN X 1 23.024 49.373 199.009 0.00 0.00 ZN

And a standard .pdb example is:
https://c4science.ch/source/lammps/browse/lammm-master-deprecated/tools/ch2lmp/example/1ac7.pdb;ef3e26cbf3b8817e379ef5a1916cdf6abee6e345
REMARK original generated coordinate pdb file
ATOM 1 P ADE 1 3.759 0.457 0.076 0.00 0.00 A
ATOM 2 O1P ADE 1 2.779 0.279 0.171 0.00 0.00 A
ATOM 3 O2P ADE 1 4.489 -0.205 0.246 0.00 0.00 A
ATOM 4 O5' ADE 1 3.969 1.972 -0.393 1.00 0.00 A
ATOM 5 C5' ADE 1 5.044 2.279 -1.290 1.00 0.00 A
ATOM 6 H5' ADE 1 4.671 2.290 -2.305 1.00 0.00 A
ATOM 7 H5'' ADE 1 5.816 1.501 -1.189 1.00 0.00 A

Do you know, or have any guesses, why atoms are still not being recognized?

xyz@hpc1:~/GBE_L1_PCA_test1$ svd GBfinal_test2.pdb z.dcd
svd: Aligning...
Subset alignment with 0 atoms converged to -nan rmsd after 1 iterations.
svd: Extracting coordinates...
svd: Allocating estimated 0.000 GB for 0 x 12 SVD
** On entry to DGESVD parameter number 6 had an illegal value
Error code from size request to dgesvd was -6

Best,
Gustavo

[agrossfield]

Like I said twice before, the default selection is 'name == "CA"' -- are there atoms called CA in your molecule? Since you're not getting anything selected, I suspect not. The atom names you're showing don't look like a protein, which means that selection wouldn't make much sense anyway. You need to figure out what atoms you want included in your PCA calculation, and write yourself a selection that gets those. For example, the selection 'backbone' will give you protein or nucleic acid backbone atoms, assuming the usual naming conventions. What you're showing looks like zinc and oxygens, so if you want just the ZN you could use 'name == "ZN"', or if you want everything you could use the selection string 'all' (though that's usually a bad idea).

You can test your selection string using the tool model-select, eg with the command
model-select --selection 'YOUR SELECTION STRING HERE' GBfinal_test2.pdb

This will print out the atoms that match your selection in a xml-ish format, and will let you figure out what atoms you're actually grabbing.

One of my previous responses has the link to the page description the LOOS selection syntax.

[gustavmf1]

Hi Alan,

Thank you very much for your last answer and advice.

I think the .pdb file has been properly read, and we are progressing.
However, after running:

user@hpc1:~/GBE_L1_PCA_test2$ svd GBF.pdb traj.dcd -A name == "ZN" -S name == "ZN"
terminate called after throwing an instance of 'std::runtime_error'
what(): Error- trajectory filename must end in an extension or the filetype must be explicitly specified
Aborted (core dumped)

We are facing this new problem: trajectory filename must end in an extension or the filetype must be explicitly specified.
Could you please give me some idea of how to overcome this error?
I've been looking into the code below to try to find a solution regarding the suffix.empty() (line 107)
https://grossfieldlab.github.io/loos/sfactories_8cpp_source.html

Thank you in advance,
Best,
Gustavo

[agrossfield]

I still see 2 problems. First, the options (-A name == "ZN" -S name == "ZN" should go before the required arguments. Second, you need to embed them in quotes, just like in the examples I've given you repeatedly. So, it should be something like -A 'name == "ZN"' -S 'name == "ZN"' or the shell won't interpret the whole selection string as a single argument. Instead, it's seeing the option as "-A name", followed by arguments "==" and "ZN", which it's trying to interpret as trajectories, and thus is failing.

Thus, the full correct command line should probably be:

svd -A 'name == "ZN"' -S 'name == "ZN"' GBF.pdb traj.dcd

I pointed this out above -- all of my examples have quotes around the arguments, as do all of the examples on the selection string web page (which also has a whole section about the need to quote them). I confess I'm getting a bit frustrated, because it feels like you're not bothering to read the things I'm sending you.

[gustavmf1]

Alan, thank you so much!

Yes, now I understood. I didn't realize about the quotes and the order, I was really focused in getting the correct .pdb file.
It ran smooth:
user@hpc1:~/GBE_L1_PCA_test2$ svd -A 'name == "ZN"' -S 'name == "ZN"' GBF.pdb traj.dcd
svd: Aligning...
Subset alignment with 4682 atoms converged to 1.23169e-12 rmsd after 1 iterations.
svd: Extracting coordinates...
svd: Allocating estimated 1.470 GB for 14046 x 1 SVD
svd: SVD requests 449474 extra space for a grand total of 1.47338 GB
svd: Calculating SVD...
svd: Done! Calculation took 0s
svd: Writing results...
svd: done!

Sorry for the many questions, I really didn't mean to be redundant.
Grateful!
Best wishes,
Gustavo