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