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 package ffx.algorithms.optimize;
39
40 import ffx.algorithms.cli.ManyBodyOptions;
41 import ffx.potential.ForceFieldEnergy;
42 import ffx.potential.MolecularAssembly;
43 import ffx.potential.bonded.AminoAcidUtils;
44 import ffx.potential.bonded.Atom;
45 import ffx.potential.bonded.Residue;
46 import ffx.potential.bonded.Rotamer;
47 import ffx.potential.bonded.RotamerLibrary;
48 import ffx.potential.openmm.OpenMMEnergy;
49 import ffx.potential.parameters.ForceField;
50 import ffx.potential.parameters.TitrationUtils;
51 import ffx.potential.parsers.PDBFilter;
52
53 import java.io.File;
54 import java.util.HashSet;
55 import java.util.List;
56 import java.util.Set;
57 import java.util.logging.Logger;
58
59 import static ffx.potential.bonded.RotamerLibrary.applyRotamer;
60
61 public class TitrationManyBody {
62
63 private static final Logger logger = Logger.getLogger(TitrationManyBody.class.getName());
64
65 private final ForceField forceField;
66 private final List<Integer> resNumberList;
67 private final double pH;
68 private final String filename;
69 private PDBFilter protonFilter;
70 private ForceFieldEnergy potentialEnergy;
71 private ManyBodyOptions manyBodyOptions;
72
73 public Set<Atom> getExcludeAtoms() {
74 return excludeAtoms;
75 }
76
77 public void setExcludeAtoms(Set<Atom> excludeAtoms) {
78 this.excludeAtoms = excludeAtoms;
79 }
80
81 private Set<Atom> excludeAtoms = new HashSet<>();
82
83 public TitrationManyBody(String filename, ForceField forceField, List<Integer> resNumberList,
84 double pH) {
85 this.filename = filename;
86 this.forceField = forceField;
87 this.resNumberList = resNumberList;
88 this.pH = pH;
89 }
90
91 public TitrationManyBody(String filename, ForceField forceField, List<Integer> resNumberList,
92 double pH, ManyBodyOptions manyBodyOptions) {
93 this.filename = filename;
94 this.forceField = forceField;
95 this.resNumberList = resNumberList;
96 this.pH = pH;
97 this.manyBodyOptions = manyBodyOptions;
98 }
99
100 public MolecularAssembly getProtonatedAssembly() {
101 MolecularAssembly protonatedAssembly = new MolecularAssembly(filename);
102 protonatedAssembly.setForceField(forceField);
103 File structureFile = new File(filename);
104 protonFilter = new PDBFilter(structureFile, protonatedAssembly, forceField,
105 forceField.getProperties(), resNumberList);
106 protonFilter.setRotamerTitration(true);
107 protonFilter.readFile();
108 protonFilter.applyAtomProperties();
109 protonatedAssembly.finalize(true, forceField);
110 potentialEnergy = ForceFieldEnergy.energyFactory(protonatedAssembly);
111 protonatedAssembly.setFile(structureFile);
112 double proteinDielectric = 1.0;
113 boolean tanhCorrection = false;
114 try {
115 proteinDielectric = forceField.getDouble("SOLUTE_DIELECTRIC", 1.0);
116 tanhCorrection = forceField.getBoolean("TANH_CORRECTION", true);
117 } catch (Exception e) {
118 logger.info("Protein Dielectric or Tanh Correction is Null");
119 }
120
121 TitrationUtils titrationUtils;
122 titrationUtils = new TitrationUtils(protonatedAssembly.getForceField(), proteinDielectric, tanhCorrection);
123 titrationUtils.setRotamerPhBias(298.15, pH);
124 for (Residue residue : protonatedAssembly.getResidueList()) {
125 String resName = residue.getName();
126 if (resNumberList.contains(residue.getResidueNumber())) {
127 if (resName.equalsIgnoreCase("ASH") || resName.equalsIgnoreCase("GLH")
128 || resName.equalsIgnoreCase("LYS") || resName.equalsIgnoreCase("HIS")
129 || resName.equalsIgnoreCase("CYS")) {
130 residue.setTitrationUtils(titrationUtils);
131 }
132 }
133 }
134 if (potentialEnergy instanceof OpenMMEnergy openMMEnergy) {
135 boolean updateBondedTerms = forceField.getBoolean("TITRATION_UPDATE_BONDED_TERMS", true);
136 openMMEnergy.getSystem().setUpdateBondedTerms(updateBondedTerms);
137 }
138 potentialEnergy.energy();
139 return protonatedAssembly;
140 }
141
142 public MolecularAssembly[] getProtonatedAssemblies() {
143 logger.info(" Adding protons that titrate to molecular assemblies");
144 MolecularAssembly molecularAssembly = getProtonatedAssembly();
145 List<Character> altLocs = protonFilter.getAltLocs();
146 int locs = 1;
147 if (altLocs != null) {
148 locs = altLocs.size();
149 for (int i = 0; i < locs; i++) {
150 if (altLocs.get(i) == null) {
151 altLocs.remove(altLocs.get(i));
152 }
153 }
154 }
155 MolecularAssembly[] molecularAssemblies = new MolecularAssembly[locs];
156 molecularAssemblies[0] = molecularAssembly;
157 for (int i = 1; i < altLocs.size(); i++) {
158 logger.info(filename);
159 MolecularAssembly newAssembly = new MolecularAssembly(filename);
160 newAssembly.setForceField(forceField);
161 File structureFile = new File(filename);
162 protonFilter = new PDBFilter(structureFile, newAssembly, forceField,
163 forceField.getProperties(), resNumberList);
164 protonFilter.setRotamerTitration(true);
165 protonFilter.setAltID(newAssembly, altLocs.get(i));
166 protonFilter.readFile();
167 protonFilter.applyAtomProperties();
168 newAssembly.finalize(true, forceField);
169 potentialEnergy = ForceFieldEnergy.energyFactory(newAssembly);
170
171 double proteinDielectric = forceField.getDouble("SOLUTE_DIELECTRIC", 1.0);
172 boolean tanhCorrection = forceField.getBoolean("TANH_CORRECTION", true);
173 TitrationUtils titrationUtils = new TitrationUtils(molecularAssembly.getForceField(),
174 proteinDielectric, tanhCorrection);
175 titrationUtils.setRotamerPhBias(298.15, pH);
176 for (Residue residue : molecularAssembly.getResidueList()) {
177 String resName = residue.getName();
178 if (resNumberList.contains(residue.getResidueNumber())) {
179 if (resName.equalsIgnoreCase("ASH") || resName.equalsIgnoreCase("GLH")
180 || resName.equalsIgnoreCase("LYS") || resName.equalsIgnoreCase("HIS")
181 || resName.equalsIgnoreCase("CYS")) {
182 residue.setTitrationUtils(titrationUtils);
183 }
184 }
185 }
186 potentialEnergy.energy();
187 molecularAssemblies[i] = newAssembly;
188 }
189 return molecularAssemblies;
190 }
191
192
193 public boolean excludeExcessAtoms(Set<Atom> excludeAtoms, int[] optimalRotamers,
194 List<Residue> residueList) {
195 boolean isTitrating = false;
196 int i = 0;
197 for (Residue residue : residueList) {
198 RotamerLibrary rotamerLibrary = manyBodyOptions.getRotamerLibrary(true);
199 residue.setRotamers(rotamerLibrary);
200 String resName = residue.getName();
201 if (resName.equalsIgnoreCase("ASH") || resName.equalsIgnoreCase("GLH")
202 || resName.equalsIgnoreCase("LYS") || resName.equalsIgnoreCase("HIS")
203 || resName.equalsIgnoreCase("CYS")) {
204 Rotamer rotamer = residue.getRotamers()[optimalRotamers[i]];
205 applyRotamer(residue, rotamer);
206 if (residue.getTitrationUtils() != null) {
207 isTitrating = true;
208 AminoAcidUtils.AminoAcid3 aa3 = rotamer.aminoAcid3;
209 residue.setName(aa3.name());
210 switch (aa3) {
211 case HID, GLU -> {
212
213 Atom HE2 = residue.getAtomByName("HE2", true);
214 excludeAtoms.add(HE2);
215 }
216 case HIE -> {
217
218 Atom HD1 = residue.getAtomByName("HD1", true);
219 excludeAtoms.add(HD1);
220 }
221 case ASP -> {
222
223 Atom HD2 = residue.getAtomByName("HD2", true);
224 excludeAtoms.add(HD2);
225 }
226 case LYD -> {
227
228 Atom HZ3 = residue.getAtomByName("HZ3", true);
229 excludeAtoms.add(HZ3);
230 }
231 case CYD -> {
232
233 Atom HG = residue.getAtomByName("HG", true);
234 excludeAtoms.add(HG);
235 }
236 default -> {
237 }
238
239 }
240 }
241 }
242 i++;
243 }
244
245 setExcludeAtoms(excludeAtoms);
246
247 return isTitrating;
248 }
249
250 public boolean excludeExcessAtoms(Set<Atom> excludeAtoms, int[] optimalRotamers,
251 MolecularAssembly molecularAssembly, List<Residue> residueList) {
252 boolean isTitrating = false;
253 double proteinDielectric = 1.0;
254 boolean tanhCorrection = false;
255 try {
256 proteinDielectric = forceField.getDouble("SOLUTE_DIELECTRIC", 1.0);
257 tanhCorrection = forceField.getBoolean("TANH_CORRECTION", true);
258 } catch (Exception e) {
259 logger.info("Protein Dielectric or Tanh Correction is Null");
260 }
261 TitrationUtils titrationUtils;
262 titrationUtils = new TitrationUtils(molecularAssembly.getForceField(), proteinDielectric, tanhCorrection);
263 titrationUtils.setRotamerPhBias(298.15, pH);
264 for (Residue residue : residueList) {
265 String resName = residue.getName();
266 if (resNumberList.contains(residue.getResidueNumber())) {
267 if (resName.equalsIgnoreCase("ASH") || resName.equalsIgnoreCase("GLH")
268 || resName.equalsIgnoreCase("LYS") || resName.equalsIgnoreCase("HIS")
269 || resName.equalsIgnoreCase("CYS")) {
270 residue.setTitrationUtils(titrationUtils);
271 }
272 }
273 }
274 isTitrating = excludeExcessAtoms(excludeAtoms, optimalRotamers, residueList);
275 return isTitrating;
276 }
277
278
279 }