Convergence Time

[Emily Sprague-Klein]

Hi,

I'm new to fullrmc and am running your THF example to test out the package on my computer.  The current calculation has been running for over 24 hours (on a 16 GB memory, 4 core/8 processors laptop), and it appears to need more time before the stochastic model converges.  If I kill the python script at any time before the calculation converges on its own, will the program automatically save the last known configuration?  

Thanks,
Emily

[Bachir Aoun]

Hi Emily,

The THF example is meant to demonstrate some serious fullrmc capabilities.  Your stochastic engine can isolate preforming angular moves, bonds moves, rotations and translations about and along symmetry axes, etc.

the simulation is purposefully long and it will converge eventually but this is not the goal of the example per say. Users are meant to learn from this example and extract code snippets from it and adapt to their use cases.

The engine saves itself upon run time and this depends on saveFrequency parameters. If you stop your engine properly (not when it's saving) you should be able to load it back and resume your simulation.

you can always export a pdb snapshot of your engine state by calling Engine.export_pdb method

hope that answers your question

regards 

 

[Emily Sprague-Klein]

Hi Bachir,

Thanks so much for your response, it was very helpful.  I ended up running the THF example successfully and have began a RMC simulation on my PDF/HEXS data set on cobalt cubane complexes in water.  Attached is my run.py file, a super cell model of the cobalt cubane (no solvent) and a comparison between the current RMC simulation output with the experimental data. Would you also be able to look at the following error message I received?  Any comments and pointers would be greatly appreciated (error message is below)

Thanks so much,
Emily

 ******** RMC Output Message **********

2020-04-02 10:25:26 - fullrmc <INFO> @0 Gen:6681600 - Tr:5832345(87.290%) - Acc:5832345(87.290%) - Rem:0(0.000%) - Err:nan

2020-04-02 10:25:26 - fullrmc <INFO> @0 Runtime saving ... DON'T INTERRUPT

2020-04-02 10:25:28 - fullrmc <INFO> @0 Runtime save is successful

2020-04-02 10:25:28 - pdbparser <INFO> All records successfully exported to 'restart.pdb'

2020-04-02 10:25:28 - fullrmc <INFO> Engine @0 finishes executing all '111360' steps in 0(days) 0:18:31

2020-04-02 10:25:29 - fullrmc <ERROR> group index must be smaller than number of atoms in pdb

2020-04-02 10:25:29 - fullrmc <ERROR> None

Traceback (most recent call last):

  File "<ipython-input-1-d70430196c51>", line 1, in <module>

    runfile('C:/Users/emily/Desktop/Co cubane RMC/run.py', wdir='C:/Users/emily/Desktop/Co cubane RMC')

  File "C:\Users\emily\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 827, in runfile

    execfile(filename, namespace)

  File "C:\Users\emily\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 110, in execfile

    exec(compile(f.read(), filename, 'exec'), namespace)

  File "C:/Users/emily/Desktop/Co cubane RMC/run.py", line 694, in <module>

    bonds_all(ENGINE)

  File "C:/Users/emily/Desktop/Co cubane RMC/run.py", line 411, in bonds_all

    ENGINE.set_groups(groups)

  File "C:\Users\emily\Anaconda3\lib\site-packages\fullrmc\Engine.py", line 2039, in set_groups

    raise Exception(LOGGER.error(err))

Exception: None

[Bachir Aoun]

Emily,

Your system will never converge because you are using the wrong constraint. 

Looking at the experimental data plot you have, you must use the PairDistributionConstraint instead of PairCorrelationConstraint.

Please check the documentation to see the difference between both and make sure your experimental data respects the constraints formula.

Also the error your are getting is because the groups your are trying to set contain atom indexes that go beyond the maximum number of atoms in your system.

hope that helps

regards

[Emily Sprague-Klein]

Hi Bachir,

Thanks for the help!  I see what you mean about using the PairDistributionConstraint as I was also wondering why the output file contained a shelf feature about 1, but that makes much more sense now.  I've also adjusted the group setting so that it does not go beyond the max number of atoms in my molecule which has solved the previous issue.  

On a somewhat related note, I'm having trouble understanding the function of ENGINE.set_groups_as_molecules() and ENGINE.set_group_selector(gs) commands.  Right now I'm getting the following error message and am unsure how to assign the right number of atoms to a group such that it's recognized by the set_move_generator:

  File "C:/Users/emily/Desktop/Co cubane RMC/run.py", line 651, in molecules

    [g.set_move_generator( MoveGeneratorCollector(collection=[TranslationGenerator(amplitude=0.2),RotationGenerator(amplitude=2)],randomize=True) ) for g in ENGINE.groups]

  File "C:/Users/emily/Desktop/Co cubane RMC/run.py", line 651, in <listcomp>

    [g.set_move_generator( MoveGeneratorCollector(collection=[TranslationGenerator(amplitude=0.2),RotationGenerator(amplitude=2)],randomize=True) ) for g in ENGINE.groups]

  File "C:\Users\emily\Anaconda3\lib\site-packages\fullrmc\Core\Group.py", line 180, in set_move_generator

    raise Exception( LOGGER.error(e) )

Thanks so much!!
Emily

  

[Bachir Aoun]

Hi Emily,

I this all the error stack you are getting ? this is probably not enough to give you decisive answer but what i can tell you is that you are trying to set a collection of move generators to groups and this collection contain a RotationGenerator. A rotation is not possible if your group has less than 2 atoms !

Rotating a single atom is meaningless and this is maybe why fullrmc didn't like it.

you can probably try this

[g.set_move_generator( MoveGeneratorCollector(collection=[TranslationGenerator(amplitude=0.2),RotationGenerator(amplitude=2)],randomize=True) ) for g in ENGINE.groups if len(g)>=2]

thanks

[Emily Sprague-Klein]

Thanks so much Bachir!  That did the trick -- I applied the condition you suggested on all instances of set_move_generator requiring rotation of the entire molecule about a symmetry axis.  

I have a few additional questions regarding fullrmc and the other scripts contained within the package (sorry for the length):  

1.  I've attached the current output of the fullrmc simulation I've been working on of the cobalt cubane complex.  It appears no bond length or angle optimizations have been performed yet, and the result is purely rotational.  I think this is partly due to my changing the naming designation when  applying the BA and the B constraint (ex.:  BA_CONSTRAINT.create_angles_by_definition( anglesDefinition={"Co4O4PF6": # O-Co-O angle in cubane core )

Nevertheless, the rotational information is still very interesting.  Previously we noticed a significant change in peak intensity at outer sphere coordination distances when comparing the powder HEXS/PDF spectrum with the static crystal structure and we had  hypothesized these changes/broadening were due to changes in ligand conformation and/or rotational disorder in the pyridine ring (see attached).  Which brings me to my next question -- from looking at the attached fullrmc output, would you be inclined to assign the two broad peak features at > 8 A to ligand rotational disorder?  Is there a way to extract the simulated curve as well as the various scattering contributions?  How do I locate what the differing colored curves represent? 

2.  In plot.py, I'm receiving the following error message and it is currently outputting one graph (attached) :

 runfile('C:/Users/emily/Desktop/Co cubane RMC/plot.py', wdir='C:/Users/emily/Desktop/Co cubane RMC')

Traceback (most recent call last):

  File "<ipython-input-51-10e0682fc2b9>", line 1, in <module>

    runfile('C:/Users/emily/Desktop/Co cubane RMC/plot.py', wdir='C:/Users/emily/Desktop/Co cubane RMC')

  File "C:\Users\emily\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 827, in runfile

    execfile(filename, namespace)

  File "C:\Users\emily\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 110, in execfile

    exec(compile(f.read(), filename, 'exec'), namespace)

  File "C:/Users/emily/Desktop/Co cubane RMC/plot.py", line 19, in <module>

    PDF.plot(show=False)

  File "C:\Users\emily\Anaconda3\lib\site-packages\fullrmc\Constraints\PairDistributionConstraints.py", line 1311, in plot

    **kwargs)

  File "C:\Users\emily\Anaconda3\lib\site-packages\fullrmc\Core\Constraint.py", line 1798, in plot

    titleParams  = titleParams)

  File "C:\Users\emily\Anaconda3\lib\site-packages\fullrmc\Core\Constraint.py", line 1647, in _plot

    ax.plot(shellCenters, val, INTRA_STYLES[intraStyleIndex], label=key, **parParams )

IndexError: list index out of range  

3.  In multiplot.py, it is not recognizing 'pcf'.  I've also tried modifying that line to read totalAx.plot(shellsCenter, output["pdf"], 'k', linewidth=3.0,  markevery=25, label="total" ) which doesn't appear to solve the issue.  The error message is below: 

 runfile('C:/Users/emily/Desktop/Co cubane RMC/multiplot.py', wdir='C:/Users/emily/Desktop/Co cubane RMC')

Traceback (most recent call last):

  File "<ipython-input-52-cef71ae81735>", line 1, in <module>

    runfile('C:/Users/emily/Desktop/Co cubane RMC/multiplot.py', wdir='C:/Users/emily/Desktop/Co cubane RMC')

  File "C:\Users\emily\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 827, in runfile

    execfile(filename, namespace)

  File "C:\Users\emily\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 110, in execfile

    exec(compile(f.read(), filename, 'exec'), namespace)

  File "C:/Users/emily/Desktop/Co cubane RMC/multiplot.py", line 56, in <module>

    totalAx.plot(PDF.shellCenters, output["pcf"], 'k', linewidth=3.0,  markevery=25, label="total" )

KeyError: 'pcf'

**** Thanks so much!!!!!  I think I'm finally getting somewhere with this.  :D

[Bachir Aoun]

Emily,

I am glad this is working for you and i am happy to see that you are able to test your hypothesis and do some novel science.

About your angles, two thing might be happening:

1. If your the groups you are moving are molecular groups which means each group contain the atoms of a molecule than when performing translations or rotations upon those groups the bonds and angles will remain untouched
2. If you feel that your bonds and angles conditions are not fulfilled and getting worst during the simulation this is probably because fullrmc can't rebuild your molecules. please read this to understand how fullrmc build molecules, it's all in the pdb file. Now, assuming that your pdb file doesn't contain all the necessary information for fullrmc to parse your molecules, you wouldn't be able to set your bonds or angles using the molecular definition but you need to do that by atom indexes. Honestly i recommend you fixing your pdb file it will make your life way easier.

Regarding your structure and rotations contribution

This is great, now you will be able to appreciate the usage of frames. You will start your simulation without rotations and you will try to go as far as you can without activating the rotations. You create another frame (e.g. you name it rotations_activated). you switch to using this frame. Careful, when switching from one frame to another one, if this other one has never been used before than it will inherit all the last used frame attributes and state. Actually this is what you want so don't create your rotation frame until you really want to start doing rotations. After switching to the rotations frame, you will activate the rotations the way you are doing it now and at the end you have 2 frames to compare and contrast the results with and without rotations. You can therefore point out that the newly created features in your structure are due to the rotations ...

Regarding the plot error 

You need to manually fix multiplot.py (i will make sure to fix that in new releases)

in your multiplot.py line 56, please change

 totalAx.plot(PDF.shellCenters, output["pcf"], 'k', linewidth=3.0,  markevery=25, label="total" )

to

 totalAx.plot(PDF.shellCenters, output["total"], 'k', linewidth=3.0,  markevery=25, label="total" )

thanks

[Emily Sprague-Klein]

Ok, great thanks!  I've made those changes to multiplot,py and it looks good.  Regarding the pdbparser, my output file does not contain a number in columns 73-76 for the sequence number -- is it okay if that cell is empty or is there a value that I can manually type in (e.g. 1 perhaps?) 

I've found your tutorial on using frames, and will rewrite parts of my run.py file to include these changes.  Thanks!

[Bachir Aoun]

If you read here you will see that fullrmc groups atoms into molecules by parsing pdb file residue name, sequence number and segment identifier entities. 

Atoms sharing the same set of those entities will be considered as part of the same molecule. I think that answers your question, adding 1 to all atoms won't make any difference it's as good as empty. 

if you have enough variance in the residue name and segment identifier to split between molecules than having a sequence number won't be a problem, if not then you need to add the sequence number variance in your pdb file to make use of fullrmc's automatic recognition of molecules.