Unable to prepare system with charmm FF

[smar966]

Dear support team,

I have been trying to prepare and run adaptive MDs with with a charmm force field. I have a fatty acid bearing two double C=C bonds. I am failing to produce a good .psf topology and run the MD. Apparently the psf file is generated, but when I run the equilibration, ACEMD is unable to read the topology (I get many warning lines such as: #WARNING: illegal value: no parameters exist for dihedral c=c cr cr cr) I am using the .par and .rtf for my ligand generated with the SwissParam webserver.

I have tried many things to fix the issue, but I always failed. Namely: 1) Prepared the ligand mol2 file using different tools/approaches 2) Prepared Amber prepi and frcmod and implemented the amber force field -> this crashed during the building of the topology, apprently with some mismatch of atoms (I was unable to find any inconsistency among the files) 3) converted the final charmm structure.psf file to amber .prmtop using parmed, and then tried to run the MD as an amber MD -> acemd failed reading the prmtop file 4) I tried the same pocedure with a similar ligand without the two double bonds, and it wor prep_Linol_2.tar.gz ked well; when I added one of the double bonds the problem was back 5) Built the topology using native tleap from ambertools, and it was fine. then tried to run acemd with this topology, but in this case the equilibration became unstable at the first step.

I am using routine HTMD scripts to prepare my system, which have worked many times before with other ligands, except for this one. I am not an expert with the charmm force field, and I don't know what else to try. Can you help me find and solve the issue?

If you want to check my files, please download them with the link: https://drive.google.com/file/d/1BQxPNeyujdzw5dUTU3WJM2G_nWWStNjt/view?usp=sharing

Thank you very much in advance! Best regards, Sergio

[stefdoerr]

Hi Sergio, sorry for the delay, I was out on holidays.

Looking at the issue I am quite sure it's because of the = in the atom type name C=C. We never encountered an atom type consisting of anything other than letters and numbers so currently the parameters for that dihedral get commented out in the final parameters file:

!Following lines added from Linol_ac.par
CO2M CR   CR   CR       0.150  3     0.00 
CO2M CR   CR   HCMM    -0.070  3     0.00 
O2CM CO2M CR   CR       0.631  2   180.00 
O2CM CO2M CR   HCMM    -0.053  3     0.00 
CR   CR   CR   CR       0.051  1     0.00 
CR   CR   CR   CR       0.341  2   180.00 
CR   CR   CR   CR       0.166  3     0.00 
CR   CR   CR   HCMM     0.320  1     0.00 
CR   CR   CR   HCMM    -0.315  2   180.00 
CR   CR   CR   HCMM     0.132  3     0.00 
!CR   CR   CR   C=C     -0.147  1     0.00 
!CR   CR   CR   C=C      0.219  2   180.00 
!CR   CR   CR   C=C      0.292  3     0.00 
!CR   CR   C=C  C=C     -0.247  1     0.00 

See exclamation mark in front of all parameters having a =.

You can quickly fix it by renaming the atom type to something more conventional.

I am not 100% sure I want to make a fix for that in HTMD until I know what symbols are acceptable atom types in CHARMM. If you have any docs showing what characters are allowed it will be very appreciated and I will make a fix for it.

[smar966]

Hi Stefan, I also just came back from vacation. Thanks for your response. We will try to solve the issue by renaming the atom types and I'll let you know. Regards

[smar966]

Hi Stefan,

I modified the name of C=C to CE1 in the par file for the ligand (CE1 is listed as an alkene carbon in the par_all36m_prot.prm file), but I still got no luck running the equilibration. This time I obtained only one error line in the log.txt file: ERROR: file src/mdio/topo.cpp line 368: Conflict in redefined bond parameters for ce1 ce1

I'm not sure how I could rename the atom type better

[smar966]

Hello,

I finally managed to have it running. I manually renamed the atom type previously described as C=C to CG2D1 and it worked. When I renamed to CE1 it still did not work. It was a struggle for me to find the correct atom type, as I am not very familiar with the CHARMM force field, and the documentation is really scarce on this topic. Please consider fixing the issue in HTMD. Thanks for your help!

[vas2201]

I do have same issue with CHARMM force field. https://github.com/Acellera/htmd/issues/1050. Have you tried using amber force field for ligand ? Vas

[smar966]

Hi. Yes, with Amber I had no problems. The issue was really the naming of the atoms that came from the SwissParam website.

On Sun, 5 Nov 2023 at 16:22, Srinivasa Penumutchu @.***> wrote:

I do have same issue with CHARMM force field.

1050 https://github.com/Acellera/htmd/issues/1050.

Have you tried using amber force field for ligand ? Vas

— Reply to this email directly, view it on GitHub https://github.com/Acellera/htmd/issues/1069#issuecomment-1793767165, or unsubscribe https://github.com/notifications/unsubscribe-auth/AH7PLASVAZNRQOY64YJFZVLYC6OMFAVCNFSM6AAAAAA3ST4ULOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTOOJTG43DOMJWGU . You are receiving this because you modified the open/close state.Message ID: @.***>

[stefdoerr]

Fixed in version 2.3.8. Sorry for the delay on this, we don't really work with or support CHARMM anymore so it's relatively low on priorities.