1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
|
//
// Copyright (C) 2024-2025 Niels Maeder and other RDKit contributors
// @@ All Rights Reserved @@
// This file is part of the RDKit.
// The contents are covered by the terms of the BSD license
// which is included in the file license.txt, found at the root
// of the RDKit source tree.
//
#include <cmath>
#include <RDGeneral/test.h>
#include <catch2/catch_all.hpp>
#include <GraphMol/SmilesParse/SmilesParse.h>
#include <GraphMol/FileParsers/FileParsers.h>
#include <GraphMol/ForceFieldHelpers/UFF/UFF.h>
#include <ForceField/MMFF/Params.h>
#include <ForceField/MMFF/BondStretch.h>
#include <GraphMol/MolTransforms/MolTransforms.h>
#include "FFConvenience.h"
using namespace RDKit;
TEST_CASE("Test empty force field") {
auto mol =
"CCCO |(-1.28533,-0.0567758,0.434662;-0.175447,0.695786,-0.299881;0.918409,-0.342619,-0.555572;1.30936,-0.801512,0.71705)|"_smiles;
REQUIRE(mol);
SECTION("basics") {
auto forceField = ForceFieldsHelper::createEmptyForceFieldForMol(*mol);
REQUIRE(forceField);
forceField->initialize();
CHECK(forceField->minimize() == 0);
CHECK(forceField->calcEnergy() == 0.0);
REQUIRE(forceField->numPoints() == mol->getNumAtoms());
REQUIRE(forceField->positions().size() == mol->getNumAtoms());
auto dist = forceField->distance(0, 1);
auto pos1 = mol->getConformer().getAtomPos(0);
auto pos2 = mol->getConformer().getAtomPos(1);
CHECK(dist == (pos2 - pos1).length());
}
SECTION("add contrib and minimize") {
auto forceField = ForceFieldsHelper::createEmptyForceFieldForMol(*mol);
forceField->initialize();
auto params = ForceFields::MMFF::MMFFBond{6.0, 100.0};
auto *contrib = new ForceFields::MMFF::BondStretchContrib(forceField.get());
contrib->addTerm(0, 1, ¶ms);
forceField->contribs().push_back(ForceFields::ContribPtr(contrib));
CHECK(forceField->minimize() == 0);
CHECK(std::round(forceField->distance(0, 1)) == 100);
auto pos1 = mol->getConformer().getAtomPos(0);
auto pos2 = mol->getConformer().getAtomPos(1);
auto dist = (pos2 - pos1).length();
CHECK(std::round(dist) == 100);
}
}
TEST_CASE("github #7901") {
#if 1
auto mb = R"CTAB(Acetone, enolate form
RDKit 3D
9 8 0 0 0 0 0 0 0 0999 V2000
-1.2143 -0.3494 0.0962 C 0 0 0 0 0 0 0 0 0 0 0 0
0.1072 0.3398 0.0041 C 0 0 0 0 0 0 0 0 0 0 0 0
0.1682 1.7287 0.0554 O 0 0 0 0 0 1 0 0 0 0 0 0
1.2210 -0.3737 -0.1263 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.2311 -1.0203 0.9811 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.3951 -0.9468 -0.8225 H 0 0 0 0 0 0 0 0 0 0 0 0
-2.0281 0.3987 0.2006 H 0 0 0 0 0 0 0 0 0 0 0 0
2.1862 0.1115 -0.1943 H 0 0 0 0 0 0 0 0 0 0 0 0
2.8861 0.1115 -0.1943 H 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0
2 3 1 0
2 4 2 0
1 5 1 0
1 6 1 0
1 7 1 0
4 8 1 0
4 9 1 0
M CHG 1 3 -1
M END)CTAB";
#else
auto mb = R"CTAB(Acetone, enolate form
RDKit 3D
9 8 0 0 0 0 0 0 0 0999 V2000
-1.2143 -0.3494 0.0962 C 0 0 0 0 0 0 0 0 0 0 0 0
0.1072 0.3398 0.0041 C 0 0 0 0 0 0 0 0 0 0 0 0
0.1682 1.7287 0.0554 O 0 0 0 0 0 1 0 0 0 0 0 0
1.2210 -0.3737 -0.1263 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.2311 -1.0203 0.9811 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.3951 -0.9468 -0.8225 H 0 0 0 0 0 0 0 0 0 0 0 0
-2.0281 0.3987 0.2006 H 0 0 0 0 0 0 0 0 0 0 0 0
2.1862 0.1115 -0.1943 H 0 0 0 0 0 0 0 0 0 0 0 0
2.8861 -1.3115 -0.1943 H 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0
2 3 1 0
2 4 2 0
1 5 1 0
1 6 1 0
1 7 1 0
4 8 1 0
4 9 1 0
M CHG 1 3 -1
M END)CTAB";
#endif
v2::FileParsers::MolFileParserParams params;
params.removeHs = false;
auto mol = v2::FileParsers::MolFromMolBlock(mb, params);
REQUIRE(mol);
auto &conf = mol->getConformer();
std::unique_ptr<ForceFields::ForceField> ff{UFF::constructForceField(*mol)};
ff->initialize();
auto needsMore = ff->minimize(200);
CHECK(!needsMore);
CHECK((conf.getAtomPos(7) - conf.getAtomPos(8)).length() > 1.80);
CHECK(MolTransforms::getAngleDeg(conf, 1, 3, 8) > 110);
CHECK(MolTransforms::getAngleDeg(conf, 1, 3, 7) > 110);
CHECK(MolTransforms::getAngleDeg(conf, 7, 3, 8) > 110);
}
|
