Angle construction error after upgrading to most recent version

[micahwelsch]

Hello,

I have a polymer system that I constructed by

1. Building the molecule using Avogadro
2. Editing the atom identifiers for later step involving topology file
3. Creating a bulk system using packmol
4. Writing the topology file from scratch
5. Using a simple .tcl script from the GOMC examples that I edited to match the residues and other identifiers I set in the PDBs

After upgrading to the most recent version of GOMC, I'm running into an issue where an angle in the furan ring of my polymer molecule that was previously okay is sending the "Cannot construct angle FC3 FC2 FO !" error at the start of a simulation.

To reproduce this error I have been using the GPU executable to run in the NPT ensemble
./GOMC_GPU_NPT in-eq.conf

I expected this NPT system to condense down to a bulk density and fluctuate at the temperature and pressure set as an amorphous polymer.

Image of angle in question, provided by Gregory Schwing on Gitter page

I'm attaching a zip file with all relevant files.

Software and System Details

• OS: Linux
• Ensemble: NPT/GEMC (I also tested a previously functioning GEMC input I was using and it's returning the same new error)
• Code version: 2.70

ptf-angle-issue.zip

Additional info/context
The SMILE code for this pentamer of PTF I'm providing as an example is as follows (I believe):

O=C(O)c5ccc(C(=O)OCCCOC(=O)c4ccc(C(=O)OCCCCOC(=O)c3ccc(C(=O)OCCCC(=O)OCc2ccc(C(=O)OCCCOC(=O)c1ccc(C(=O)OCCCO)o1)o2)o3)o4)o5

[LSchwiebert]

Micah, thank you for your bug report. Can you please answer a few questions? First, I assume from your report that this was working correctly on an earlier version of the code. If so, for which version of GOMC was this working? Also, have you tried using the latest development branch? It has many bug fixes and, to the best of my knowledge, does not introduce any new bugs. I don't know that any of the changes will address your problem.

[micahwelsch]

Hi! I believe I was on version 2.4. I haven't tried the latest development branch. I can give it a shot.

[micahwelsch]

Update, the dev branch still yields the same error. From what I understand from previous commits, it could be related to the CBMC trial moves or a limitation with how GOMC handles the angles of planar molecules (i.e. the angles I have set in my forcefield file)? Could either of those be potential causes for the problem?

[LSchwiebert]

Thanks for the update. That's very helpful. It could be that the problem is caused by something you mentioned. We will look into it asap and let you know.

[jpotoff]

@GregorySchwing @LSchwiebert I suggest we take a look at the PSF parser, since we were recently doing a lot of work on how we were identifying rings in the system to address performance issues.

[micahwelsch]

I've gone ahead and tried building another polymer I'm interested in as a test and I'm not sure if it helps but it seems as if the issue is specifically being caused by angles in a ring involving an oxygen in that ring.

For example this dimer of PDD is getting a "Cannot construct angle" error for the O-C-O and C-C-O angles in the rings.
I'm attaching a screenshot of the structure with those angles highlighted.

[msoroush]

@micahwelsch I looked at your system and it seems that you angles are very flexible (angle force constants are very small). This would cause some issue with GOMC CD-CBMC algorithm. For any move that uses CD-CBMC, such as swap, intra-swap, regrowth, etc, when we sample angles in branching point, we fix the angle that belongs to ring (e.g. in the above picture O3-C6-O4) and we let crancshaft move samples this angle. This would cause issue for very flexible angle, because we cannot enforce the angle in the ring, due to geometrical constrain.
For example, in one of the branching point in your system (EsC-FC2-FO), the angle parameters are:

ANGLES
!
!V(angle) = Ktheta(Theta - Theta0)**2
! Ktheta = K * Kb
EsC  FC2  FO        70           120
FC3  FC2  EsC       85          136.7
FC3  FC2  FO        70           110.6

If we sample EsC FC2 FO to be 130 and FC3 FC2 EsC to be 140, geometry dictates that maximum angle for FC3 FC2 FO must be 90, but the actual angle in the ring that we want to keep fix is around 98. So, that's why we see the GOMC error.

If you modify the equilibrium angles that involves an atom in the ring (X-FC2-X in your system), so that the sum of them is close to 360, you should be fine.

ANGLES
!
!V(angle) = Ktheta(Theta - Theta0)**2
! Ktheta = K * Kb
EsC  FC2  FO        70           120
FC3  FC2  EsC       85          129.4
CaC  FC2  FO        70          120
FC3  FC2  CaC       85          129.4
FC3  FC2  FO        70           110.6

[micahwelsch]

Interesting. I will give this a shot. To double check probably a silly question, the Ktheta that I provide should be in kcal/mol correct? Or do I provide it in kelvin and the software converts that automatically? What I've been using and what is in the example above is in kcal/mol.

[msoroush]

Interesting. I will give this a shot. To double check probably a silly question, the Ktheta that I provide should be in kcal/mol correct? Or do I provide it in kelvin and the software converts that automatically? What I've been using and what is in the example above is in kcal/mol.

I believe you are using Mie parameter file, which force units are in kelvin. That's why angels are very flexible.
Here you can find more information regarding the difference in parameter files and how you can set it in control file.

[micahwelsch]

I fixed the energy units and still found it necessary to adjust the angles as you suggested @msoroush to get the sum closer to 360. For a polymer system with a single unit long chain (monomer) I'm getting expected results.

But with the 5 unit long chain (5 monomers) that I used as my initial example in this thread I'm still getting the angle construction error.

The system can't seem to compress below a particular density and gives the angle construction error when it does manage to. I'm not entirely sure why that might be. Here's the inputs that are giving this result.
continued-angle-const-issue.zip

[msoroush]

1. @micahwelsch I am working on a fix for your issue, and it's under the review. If you want, you can give it Fix to the issue #338 and issue #335 #339 a try.
2. I notice that you are using 10 A cutoff for electrostatic calculation (Ewald summation method). If you increase rCutCoulomb you can improve GOMC performance. Please take a look at discussion in issue Long time - 1,000 MC cycles GOMC vs Towhee - binary mixture #323
3. Regarding the density issue in your calculation, the 5 units long polymer would take longer to equilibrate compare to 1 unit polymer. Specially for atoms in middle ring. In general, longer the molecule becomes, it will be harder for MC to sample molecular conformation of groups in the middle. Moves like Regrowth and crankshaft can help to sample molecular conformation, but these moves also have their own limitations.

Please let me know if #339 patch fixes your issue.