Skip to content

Commit

Permalink
Implemented code to support the Alexandria force field.
Browse files Browse the repository at this point in the history
The Alexandria Chemistry Toolkit for building force fields
will be released shortly. It relies on the OpenBabel software to
support it for generating atom types. In this patch the
forcefieldalexandria is introduced as a child class to forcefieldgaff
in an effort to change as little code as possible.

There are many changes to bondtyp.txt which admittably will be hard to
check. These will affect any force field typing in OpenBabel.
  • Loading branch information
dspoel committed Jan 16, 2024
1 parent 32cf131 commit f186a81
Show file tree
Hide file tree
Showing 10 changed files with 816 additions and 15 deletions.
1 change: 1 addition & 0 deletions data/CMakeLists.txt
@@ -1,6 +1,7 @@
# Open Babel data files

set(to_install
alexandria.prm
atomization-energies.txt
atomtyp.txt
babel_povray3.inc
Expand Down
509 changes: 509 additions & 0 deletions data/alexandria.prm

Large diffs are not rendered by default.

160 changes: 151 additions & 9 deletions data/bondtyp.txt
Expand Up @@ -31,49 +31,178 @@
#0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
[X2,X3]1[#6]([#7]2)[#6][#6][#6]2[X2,X3][#6]([#7]3)[#6][#6][#6]3[X2,X3][#6]([#7]4)[#6][#6][#6]4[X2,X3][#6]([#7]5)[#6][#6][#6]51 0 1 2 1 2 1 1 3 1 3 4 2 4 5 1 5 2 1 5 6 2 6 7 1 7 8 2 7 9 1 9 10 2 10 11 1 11 8 1 11 12 2 12 13 1 13 14 1 13 15 2 15 16 1 16 17 2 17 14 1 17 18 1 18 19 2 19 20 1 19 21 1 21 22 2 22 23 1 23 20 2

# N-methylmethanesulfinimidic-acid
[#6X4][#7X2][#16X3] 0 1 1 1 2 2
# amino--diaminomethylideneamino-methylidene--ethylazanium and analogs
[#6X3]([#7X3])([#7X3])([#7X2][#6X3]([#7X3])([#7X3]([#6]))) 0 1 1 0 2 1 0 3 2 3 4 1 4 5 1 4 6 2 6 7 1
# cyanoimino-methylene-azanide N#C-N-C#N
[#7X1][#6X2][#7X2][#6X2][#7X1] 0 1 3 1 2 1 2 3 1 3 4 3
# fulvene analogs
[#6X3]1([#6X3]([#1])([#1]))([#6X3][#6X3][#6X3][#6X3]1) 0 1 2 1 2 1 1 3 1 0 4 1 4 5 2 5 6 1 6 7 2 0 7 1
[#6X3]1([#7X2])([#6X3][#6X3][#6X3][#6X3]1) 0 1 2 0 2 1 2 3 2 3 4 1 4 5 2 0 5 1
# CR3-N=N-CH=CH
[#6X4][#7X2][#7X2][#6X3][#6X3] 0 1 1 1 2 2 2 3 1 3 4 2
#CR2=N-N=CH2
[#6X3][#7X2][#7X2][#6X3]([#1])([#1]) 0 1 2 1 2 1 2 3 2 3 4 1 3 5 1
# C2=N-NO2
[#6X3][#7X2][#7X3]([#8X1])([#8X1]) 0 1 2 1 2 1 2 3 1 2 4 2
#1-3-oxazole analogs
[#8X2]1[#6X3][#7X2][#6X3][#6X3]1 0 1 1 1 2 2 2 3 1 3 4 2 0 4 1
#N=C=N-C#N
[#7X1][#6X2][#7X2][#6X2][#7X1] 0 1 2 1 2 2 2 3 1 3 4 3
#pyridinium analogs
[#7X3]1([#6])([#6X3][#6X3][#6X3][#6X3][#6X3]1) 0 1 1 0 2 1 2 3 2 3 4 1 4 5 2 5 6 1 0 6 2
[#7X3]1([#1])([#6X3][#6X3][#6X3][#6X3][#6X3]1) 0 1 1 0 2 2 2 3 1 3 4 2 4 5 1 5 6 2 0 6 1
#Pyrrol-imine analogs
[#6X3]1([#7X2][#1,#6])([#7X2][#6X3][#6X3][#6X3]1) 0 1 2 1 2 1 0 3 1 3 4 2 4 5 1 5 6 2 0 6 1
#1H-pyrimidin-2-one analogs
[#6X3]1([#8X1])([#7X3][#6X3][#6X3][#6X3][#7X2]1) 0 1 2 0 2 1 2 3 1 3 4 2 4 5 1 5 6 2 0 6 1
# 5-methylideneimidazol-2-one analogs, 5-methylimidazol-4-one analogs
[#6X3]1([#8X1])([#7X3][#6X3]([#6X3])[#6X3][#7X2]1) 0 1 2 0 2 1 2 3 1 3 4 2 3 5 1 5 6 2 0 6 1
[#6X3]1([#8X1])([#7X2][#6X3][#7X2][#6X3]1) 0 1 2 0 2 1 2 3 2 3 4 1 4 5 2 0 5 1
# C/N-N=CH2
[#6X3]([#1])([#1])[#7X2][#7X2,#7X3,#6X3,#6X4] 0 1 1 0 2 1 0 3 2 3 4 1
# Imines
[#6X3]([#7X2][#1])([#7X3][#6X3]([#1])[#6X3]) 0 1 2 1 2 1 0 3 1 3 4 1 4 5 1 4 6 2
[#6X3]([#7X2][#1])([#7X2][#6X3]([#1])[#6X3]) 0 1 2 1 2 1 0 3 1 3 4 2 4 5 1 4 6 1
[#6X3]([#7X2][#1,#6])([#7X3])([*]) 0 1 2 1 2 1 0 3 1 0 4 1
# Azide
[#7D2][#7D2^1][#7D1] 0 1 2 1 2 2
# Nitro
[#8D1][#7D3^2]([#8D1])* 0 1 2 1 2 2 1 3 1

# [#8D1][#7D3^2]([#8D1])* 0 1 2 1 2 2 1 3 1; old code
[#8D1][#7D3^2]([#8D1])* 0 1 1 1 2 2 1 3 1
# nitronium
[#8X1][#7X2]([#8X1]) 0 1 2 1 2 2
# Dimethyl(methylidene)azanium and analogs
[#7X3]([#6X4])([#6X4])([#6X3]([#1])([#1])) 0 1 1 0 2 1 0 3 2 3 4 1 3 5 1
[#7X3]([#6X3]([#7X3])([#7X3]))([#6X4])([#6X4]) 0 1 2 1 2 1 1 3 1 0 4 1 0 5 1
# N2O
[#7X1][#7X2][#8X1] 0 1 3 1 2 1
# nitroso R-N=O
[#8X1][#7X2][#7X2][#6X3] 0 1 2 1 2 1 2 3 2
[#8X1][#7X2]([*;!#6]) 0 1 2 1 2 1
[#8X1][#7X2]([#6X3]([*])([*])) 0 1 2 1 2 1 2 3 1 2 4 1
# H-N=P(-O)3
[#7X2]([#1])([#15X4]([#8X1])([#8X1])([#8X1])) 0 1 1 0 2 2 2 3 1 2 4 1 2 5 1
# NR2-C=C-C=O
[#7X3][#6X3][#6X3][#6X3][#8X1] 0 1 1 1 2 2 2 3 1 3 4 2
# NR2-C=C-C
[#7X3][#6X3][#6X3]([#6X4,#1])[#6X2] 0 1 1 1 2 2 2 3 1 2 4 1
# C=CH-NR-CR=O
[#6X3][#6X3]([#1])[#7X3][#6X3][#8X1] 0 1 2 1 2 1 1 3 1 3 4 1 4 5 2
# bis-methylamino-methylidene-methylazaniumyl--sulfate analogs
[#7X3]([#8X2])([#6X4])([#6X3]([#7X3])([#7X3])) 0 1 1 0 2 1 0 3 2 3 4 1 3 5 1
# NO2H=CH-CR3
[#7X3]([#8X2])([#8X1])([#6X3][#6X4]) 0 1 1 0 2 1 0 3 2 3 4 1
# CR2=N-N(CR3)2
[#6X3][#7X2][#7X3]([#6X4])([#6X4]) 0 1 2 1 2 1 2 3 1 2 4 1
#1-2-4-thiadiazole
[#6X3]1([#7X2][#6X3]([#1])[#7X2][#16X2]1) 0 1 2 1 2 1 2 3 1 2 4 2 4 5 1 0 5 1
# cyanoimino-methylene-azanide N#C-N-C#N
#[#7X1][#6X2][#7X2][#6X2][#7X1] 0 1 3 1 2 1 2 3 1 3 4 3

# CR2=CH2
[#6X3]([#1])([#1])[#6X3] 0 1 1 0 2 1 0 3 2

#R-C#N-S
[#16X1][#7X2][#6X2] 0 1 1 1 2 3
# R(RO)S=O
[#16X3]([#8X1])([#8X2])([*!#8]) 0 1 2 0 2 1 0 3 1
#R-N=S=O
[#8X1][#16X2][#7X2][*] 0 1 2 1 2 2 2 3 1
[#8X1][#16X4]([#7X2][#1])([*])([*]) 0 1 2 1 2 2 2 3 1 1 4 1 1 5 1
# Sulfones
[#16D4]([#8D1])([#8D1])([*!#8])([*!#8]) 0 1 2 0 2 2 0 3 1 0 4 1
# Sulfates
[#16D4]([#8D1])([#8D1])([#8-,#8D1])([#8-,#8D1]) 0 1 2 0 2 2 0 3 1 0 4 1
# Thiosulfates
[#16D4]([#16D1])([#8D1])([#8-,#8])([#8-,#8]) 0 1 2 0 2 2 0 3 1 0 4 1
#[#16D4]([#16D1])([#8D1])([#8-,#8])([#8-,#8]) 0 1 2 0 2 2 0 3 1 0 4 1
[#16D4]([#16D1])([#8D1])([#8])([#8-,#8]) 0 1 1 0 2 2 0 3 2 0 4 1
# Sulfonyl analogs
[#16X3]([#8X1])([#8X1])([#6^2,#7^2]) 0 1 2 0 2 2 0 3 2
# Sulfinic acid and analogs
[#16X3]([#8X1])([#8X1,#8X2])([#6]) 0 1 2 0 2 1 0 3 1
# Sulfoxides
[#16D3]([#8D1])([*!#8])([*!#8]) 0 1 2 0 2 1 0 3 1
# Sulfite
[#16D3]([#8D1])([#8D1-])([#8D1-]) 0 1 2 0 2 1 0 3 1
# Sulfur trioxide
[#16D3^2]([#8D1])([#8D1])([#8D1]) 0 1 2 0 2 2 0 3 2
[#16D3^2]([#8D1])([#8D1])([#8D1]) 0 1 2 0 2 2 0 3 2
# Sulfites
[#16D3]([#8D1])([#8])([#8]) 0 1 2 0 2 1 0 3 1
# 1-Sulfinylethane and analogs C=S=O
[#8X1][#16X2][#6X3] 0 1 2 1 2 2
#H2C=S=O
[#16X4]([#8X1])([#6X3]([#1])([#1])) 0 1 2 0 2 2 2 3 1 2 4 1
# Sulfonic acid and sulfonic ester
[#16X4]([#8X1])([#8X1])([#8])* 0 1 2 0 2 2 0 3 1 0 4 1
# Disulfur monoxide
[#16D2]([#8D1])([#16D1]) 0 1 2 0 2 2
# Sulfmonoxides
[#16D2]([#8D1])([*!#8]) 0 1 2 0 2 1
# S=N-H
[#16X2,#16X4][#7X2][#1] 0 1 2 1 2 1
# Sulfur dioxide
[#16D2]([#8D1])([#8D1]) 0 1 2 0 2 2

# S=S=S, S=S
[#16X1][#16X2][#16X1] 0 1 2 1 2 2
[#16X1][#16X1] 0 1 2
# S=CH2
[#16X2][#6X3]([#1])([#1]) 0 1 2 1 2 1 1 3 1
# dimethyl-oxido-sulfanylidene-lambda5-phosphane and analogs
[#16X1][#15] 0 1 2
# dihydroxy-sulfanylidene--lambda4-sulfane and analogs
[#16X1;!-][#16X3,#16X4] 0 1 2
# thiopyrylium
[#16X2]1[#6X3][#6X3][#6X3]([#1])[#6X3][#6X3]1 0 1 1 1 2 2 2 3 1 3 4 1 3 5 2 5 6 1 0 6 2
# dithiol-1-ium analogs
[#16X2]1[#6X3][#6X3][#6X3][#16X2]1 0 1 2 1 2 1 2 3 2 3 4 1 0 4 1
# O=S=CR2
[#8X1][#16X2][#6X3] 0 1 2 1 2 2

# 4-methylidenetriphosphole
[#6X3]1([#6X3])([#15X2][#15X2][#15X2][#6X3]1) 0 1 2 0 2 1 2 3 2 3 4 1 4 5 2 0 5 1
# C=P-H
[#6X3]1([#15X2][#1])([#6X3][#15X2][#6X3][#15X2]1) 0 1 2 1 2 1 0 3 1 3 4 2 4 5 1 5 6 2 0 6 1
[#6X3]1[#15X2][#6X3][#6X3][#6X3]1 0 1 1 1 2 2 2 3 1 3 4 2 0 4 1
[#6X3,#6X2,#16X2][#15X2][#1] 0 1 2 1 2 1
# 3-4-dihydrophosphol-2-one analogs
[#6X3]1[#15X2][#6X3][#6X4][#6X4]1 0 1 1 1 2 2 2 3 1 3 4 1 0 4 1
# C=P-C
[#6X3][#15X2][#6;!r6] 0 1 2 1 2 1
#Phosphite
[#15D3]([#8D1])([#8D1])([#8D2]) 0 1 2 0 2 2 0 3 1

#Methylphosphinate and analogs
[#15X3^2]([#8X1;!-])([#8X1;!-])([#6]) 0 1 2 0 2 2 0 3 1
#oxophosphane
#[#15D2]([#8D1])([#1]) 0 1 2 0 2 1
[#15D2]([#8D1])([#1]) 0 1 2 0 2 1
#phosphaalkene and analogs
[#15X2;!R]([#6X3,#6X4])([#6X3]) 0 1 1 0 2 2
[#15X2]([#6X3;!H,#6X4])([#6X3;!H]) 0 1 1 0 2 2
# P=O
[#15X2][#8X1] 0 1 2

#Nitrosyl Hydride
[#7D2]([#8D1])([#1]) 0 1 2 0 2 1
#Nitrite
[#8X1][#7X2][#8X2][*] 0 1 2 1 2 1 2 3 1

# chloro-hydroxy-dioxidophosphanium; not working right now
#[#15X4]([#8X1])([#8X1])([#8X2])([*]) 0 1 1 0 2 1 0 3 1 0 4 1
#[#15X4]([#8X1])([#8X1])([#17])([*]) 0 1 1 0 2 1 0 3 1 0 4 1
# Phosphone
[#15D4]([#8D1])(*)(*)(*) 0 1 2 0 2 1 0 3 1 0 4 1

[#15X4]([#8X1])([*])([*])([*]) 0 1 2 0 2 1 0 3 1 0 4 1
# diphospholes and analogs
[#6][#15X2][#15X2][#6] 0 1 1 1 2 2 2 3 1
# C-P=O
[#15X2]([#8X1])([#6]) 0 1 2 0 2 1

# Carboxylic Acid, ester, etc.
[#6D3^2]([#8D1])([#8])* 0 1 2 0 2 1 0 3 1
# Carbon dioxide
[#8D1][#6D2^1][#8D1] 0 1 2 1 2 2
# Carbon monoxide
[#8X1][#6X1] 0 1 3
# Amide C(=O)N - no negative charge on O (2aix_u1k.sdf)
[#6D3^2]([#8D1;!-])([#7])* 0 1 2 0 2 1 0 3 1
# Seleninic acid Se(=O)OH
Expand All @@ -86,14 +215,26 @@
# avoid aromatics (pdb_ligands_sdf/1yry_msg.sdf)
[CD3^2]([#16D1])([N])* 0 1 2 0 2 1 0 3 1

# C=C=C=CH2
[#6X3][#6X2][#6X2][#6X3]([#1])([#1]) 0 1 2 1 2 2 2 3 2 3 4 1 3 5 1
# allene C=C=C
# (this is problematic -- need to make sure the center carbon is sp)
[#6^2][#6D2^1][#6^2] 0 1 2 1 2 2

# ethenolate
[#6X3]([#1])([#1])[#6X3]([#8])([#1]) 0 1 1 0 2 1 0 3 2 3 4 1 3 5 1
# acetylenoate
[#6X2]([#1])[#6X2]([#8]) 0 1 1 0 2 2 2 3 1

# ene-one C=C=O
[#6^2][#6D2^1][#8D1] 0 1 2 1 2 2

# isonitrile / isocyano
[#6D1][#7D2^1]* 0 1 3 1 2 1
# NR2 in ring with hybridized carbon neighbors, do not apply to aromatic N
# valence insted of degree used to fix pdb_ligands_sdf/3dcv_55e.sdf
[Nv2R][#6v3^2][#8v2] 0 1 2 1 2 1
[Nv2R][#6v3^2][Nv2] 0 1 2 1 2 1

# if three N are present in R-N-guanidine-ish, prefer double bond to the
# non-terminal N (i.e. D2 if present)
Expand All @@ -111,3 +252,4 @@

### other potential functional groups that may (or may not) be useful to add
# imidines ( N=C/N\C=N

3 changes: 3 additions & 0 deletions data/ringtyp.txt
Expand Up @@ -12,6 +12,9 @@
# #
##############################################################################

# 3 membered rings
RINGTYP thiirene S1C=C1

#
# 5 membered rings
#
Expand Down
1 change: 1 addition & 0 deletions src/CMakeLists.txt
Expand Up @@ -45,6 +45,7 @@ set(fingerprints
)

set(forcefields
forcefields/forcefieldalexandria.cpp
forcefields/forcefieldgaff.cpp
forcefields/forcefieldghemical.cpp
forcefields/forcefieldmmff94.cpp
Expand Down
48 changes: 46 additions & 2 deletions src/bondtyper.cpp
Expand Up @@ -115,8 +115,8 @@ namespace OpenBabel

OBSmartsPattern *currentPattern;
OBBond *b1, *b2;
OBAtom *a1,*a2, *a3;
double angle, dist1, dist2;
OBAtom *a1,*a2, *a3, *a4;
double angle, angle1, angle2, dist1, dist2;
vector<int> assignments;
vector<vector<int> > mlist;
vector<vector<int> >::iterator matches, l;
Expand Down Expand Up @@ -207,6 +207,50 @@ namespace OpenBabel
}
} // thione

// Imidazole [nD3]1c=[nD3]c=c1
OBSmartsPattern imidazolium; imidazolium.Init("[#7D3][#6][#7D3][#6][#6]");

if (imidazolium.Match(mol))
{
mlist = imidazolium.GetUMapList();
for (l = mlist.begin(); l != mlist.end(); ++l)
{
a1 = mol.GetAtom((*l)[1]);
a2 = mol.GetAtom((*l)[2]);
a2->SetFormalCharge(1);

angle1 = a2->AverageBondAngle();
dist1 = a1->GetDistance(a2);

if (angle1 > 115 && angle < 150 && dist1 < 1.72) {

if ( !a1->HasDoubleBond() ) {// no double bond already assigned
b1 = a1->GetBond(a2);

if (!b1 ) continue;
b1->SetBondOrder(2);
}
}

a3 = mol.GetAtom((*l)[3]);
a4 = mol.GetAtom((*l)[4]);

angle2 = a3->AverageBondAngle();
dist2 = a3->GetDistance(a4);

// imidazolium geometries ?
if (angle2 > 115 && angle2 < 150 && dist2 < 1.72) {

if ( !a3->HasDoubleBond() ) {// no double bond already assigned
b2 = a3->GetBond(a4);

if (!b2 ) continue;
b2->SetBondOrder(2);
}
}
}
} // imidazolium

// Isocyanate N=C=O or Isothiocyanate
bool dist1OK;
OBSmartsPattern isocyanate; isocyanate.Init("[#8,#16;D1][#6D2][#7D2]");
Expand Down
39 changes: 39 additions & 0 deletions src/forcefields/forcefieldalexandria.cpp
@@ -0,0 +1,39 @@
/**********************************************************************
forcefieldalexandria.cpp - Alexandria force field.
Copyright (C) 2009 by Frank Peters <e.a.j.f.peters@tue.nl>
Copyright (C) 2006-2007 by Tim Vandermeersch <tim.vandermeersch@gmail.com>
Copyright (C) 2021 by David van der Spoel <david.vanderspoel@icm.uu.se>
This file is part of the Open Babel project.
For more information, see <http://openbabel.org/>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation version 2 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
***********************************************************************/

#include <openbabel/babelconfig.h>

#include "forcefieldalexandria.h"

#include <cstdlib>

using namespace std;

namespace OpenBabel
{
//***********************************************
//Make a global instance
OBForceFieldAlexandria theForceFieldAlexandria("Alexandria", true);
//***********************************************

} // end namespace OpenBabel

//! \file forcefieldalexandria.cpp
//! \brief Alexandria force field
50 changes: 50 additions & 0 deletions src/forcefields/forcefieldalexandria.h
@@ -0,0 +1,50 @@
/**********************************************************************
forcefieldalexandriaf.h - Alexandria force field.
Copyright (C) 2009 by Frank Peters <e.a.j.f.peters@tue.nl>
Copyright (C) 2006 by Tim Vandermeersch <tim.vandermeersch@gmail.com>
Copyright (C) 2021 by David van der Spoel <david.vanderspoel@icm.uu.se>
This file is part of the Open Babel project.
For more information, see <http://openbabel.org/>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation version 2 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
***********************************************************************/

#include <vector>
#include <string>
#include <map>

#include <openbabel/forcefield.h>
#include <openbabel/base.h>
#include <openbabel/mol.h>

#include "forcefieldgaff.h"

namespace OpenBabel
{

// Class OBForceFieldAlexandria
// class introduction in forcefieldalexandria.cpp
class OBForceFieldAlexandria: public OBForceFieldGaff
{
public:
//! Constructor
explicit OBForceFieldAlexandria(const char* ID, bool IsDefault=true) : OBForceFieldGaff(ID, IsDefault)
{
SetPrmFile("alexandria.prm");
}

}; // class OBForceFieldAlexandria

}// namespace OpenBabel

//! \file forcefieldalexandria.h
//! \brief Alexandria force field

0 comments on commit f186a81

Please sign in to comment.