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("Getting protonated 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 = 0; i < altLocs.size(); i++) {
158 if (i != 0) {
159 logger.info(filename);
160 MolecularAssembly newAssembly = new MolecularAssembly(filename);
161 newAssembly.setForceField(forceField);
162 File structureFile = new File(filename);
163 protonFilter = new PDBFilter(structureFile, newAssembly, forceField,
164 forceField.getProperties(), resNumberList);
165 protonFilter.setRotamerTitration(true);
166 protonFilter.setAltID(newAssembly, altLocs.get(i));
167 protonFilter.readFile();
168 protonFilter.applyAtomProperties();
169 newAssembly.finalize(true, forceField);
170 potentialEnergy = ForceFieldEnergy.energyFactory(newAssembly);
171 double proteinDielectric = 1.0;
172 boolean tanhCorrection = false;
173 try {
174 proteinDielectric = forceField.getDouble("SOLUTE_DIELECTRIC", 1.0);
175 tanhCorrection = forceField.getBoolean("TANH_CORRECTION", true);
176 } catch (Exception e) {
177 logger.info("Protein Dielectric or Tanh Correction is Null");
178 }
179
180 TitrationUtils titrationUtils;
181 titrationUtils = new TitrationUtils(molecularAssembly.getForceField(), proteinDielectric,tanhCorrection);
182 titrationUtils.setRotamerPhBias(298.15, pH);
183 for (Residue residue : molecularAssembly.getResidueList()) {
184 String resName = residue.getName();
185 if (resNumberList.contains(residue.getResidueNumber())) {
186 if (resName.equalsIgnoreCase("ASH") || resName.equalsIgnoreCase("GLH")
187 || resName.equalsIgnoreCase("LYS") || resName.equalsIgnoreCase("HIS")
188 || resName.equalsIgnoreCase("CYS")) {
189 residue.setTitrationUtils(titrationUtils);
190 }
191 }
192 }
193 potentialEnergy.energy();
194 molecularAssemblies[i] = newAssembly;
195 }
196 }
197 return molecularAssemblies;
198 }
199
200
201 public boolean excludeExcessAtoms(Set<Atom> excludeAtoms, int[] optimalRotamers,
202 List<Residue> residueList) {
203 boolean isTitrating = false;
204 int i = 0;
205 for (Residue residue : residueList) {
206 RotamerLibrary rotamerLibrary = manyBodyOptions.getRotamerLibrary(true);
207 residue.setRotamers(rotamerLibrary);
208 String resName = residue.getName();
209 if(resName.equalsIgnoreCase("ASH") || resName.equalsIgnoreCase("GLH")
210 || resName.equalsIgnoreCase("LYS") || resName.equalsIgnoreCase("HIS")
211 || resName.equalsIgnoreCase("CYS")){
212 Rotamer rotamer = residue.getRotamers()[optimalRotamers[i]];
213 applyRotamer(residue, rotamer);
214 if (residue.getTitrationUtils() != null) {
215 isTitrating = true;
216 AminoAcidUtils.AminoAcid3 aa3 = rotamer.aminoAcid3;
217 residue.setName(aa3.name());
218 switch (aa3) {
219 case HID, GLU -> {
220
221 Atom HE2 = residue.getAtomByName("HE2", true);
222 excludeAtoms.add(HE2);
223 }
224 case HIE -> {
225
226 Atom HD1 = residue.getAtomByName("HD1", true);
227 excludeAtoms.add(HD1);
228 }
229 case ASP -> {
230
231 Atom HD2 = residue.getAtomByName("HD2", true);
232 excludeAtoms.add(HD2);
233 }
234 case LYD -> {
235
236 Atom HZ3 = residue.getAtomByName("HZ3", true);
237 excludeAtoms.add(HZ3);
238 }
239 case CYD -> {
240
241 Atom HG = residue.getAtomByName("HG", true);
242 excludeAtoms.add(HG);
243 }
244 default -> {
245 }
246
247 }
248 }
249 }
250 i++;
251 }
252
253 setExcludeAtoms(excludeAtoms);
254
255 return isTitrating;
256 }
257
258 public boolean excludeExcessAtoms(Set<Atom> excludeAtoms, int[] optimalRotamers,
259 MolecularAssembly molecularAssembly, List<Residue> residueList){
260 boolean isTitrating = false;
261 double proteinDielectric = 1.0;
262 boolean tanhCorrection = false;
263 try {
264 proteinDielectric = forceField.getDouble("SOLUTE_DIELECTRIC", 1.0);
265 tanhCorrection = forceField.getBoolean("TANH_CORRECTION", true);
266 } catch (Exception e) {
267 logger.info("Protein Dielectric or Tanh Correction is Null");
268 }
269 TitrationUtils titrationUtils;
270 titrationUtils = new TitrationUtils(molecularAssembly.getForceField(), proteinDielectric,tanhCorrection);
271 titrationUtils.setRotamerPhBias(298.15, pH);
272 for (Residue residue : residueList) {
273 String resName = residue.getName();
274 if (resNumberList.contains(residue.getResidueNumber())) {
275 if (resName.equalsIgnoreCase("ASH") || resName.equalsIgnoreCase("GLH")
276 || resName.equalsIgnoreCase("LYS") || resName.equalsIgnoreCase("HIS")
277 || resName.equalsIgnoreCase("CYS")) {
278 residue.setTitrationUtils(titrationUtils);
279 }
280 }
281 }
282 isTitrating = excludeExcessAtoms(excludeAtoms,optimalRotamers, residueList);
283 return isTitrating;
284 }
285
286
287 }