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.potential.ForceFieldEnergy;
43 import ffx.potential.openmm.OpenMMEnergy;
44 import ffx.potential.MolecularAssembly;
45 import ffx.potential.bonded.AminoAcidUtils;
46 import ffx.potential.bonded.Atom;
47 import ffx.potential.bonded.Residue;
48 import ffx.potential.bonded.Rotamer;
49 import ffx.potential.openmm.OpenMMEnergy;
50 import ffx.potential.parameters.ForceField;
51 import ffx.potential.parameters.TitrationUtils;
52 import ffx.potential.parsers.PDBFilter;
53
54 import java.io.File;
55 import java.util.List;
56 import java.util.Set;
57 import java.util.logging.Logger;
58 import java.lang.String;
59
60 public class TitrationManyBody {
61
62 private static final Logger logger = Logger.getLogger(TitrationManyBody.class.getName());
63
64 private final ForceField forceField;
65 private final List<Integer> resNumberList;
66 private final double pH;
67 private final String filename;
68 private PDBFilter protonFilter;
69 private ForceFieldEnergy potentialEnergy;
70
71 public TitrationManyBody(String filename, ForceField forceField, List<Integer> resNumberList,
72 double pH) {
73 this.filename = filename;
74 this.forceField = forceField;
75 this.resNumberList = resNumberList;
76 this.pH = pH;
77 }
78
79 public MolecularAssembly getProtonatedAssembly() {
80 MolecularAssembly protonatedAssembly = new MolecularAssembly(filename);
81 protonatedAssembly.setForceField(forceField);
82 File structureFile = new File(filename);
83 protonFilter = new PDBFilter(structureFile, protonatedAssembly, forceField,
84 forceField.getProperties(), resNumberList);
85 protonFilter.setRotamerTitration(true);
86 protonFilter.readFile();
87 protonFilter.applyAtomProperties();
88 protonatedAssembly.finalize(true, forceField);
89 potentialEnergy = ForceFieldEnergy.energyFactory(protonatedAssembly);
90 protonatedAssembly.setFile(structureFile);
91
92 TitrationUtils titrationUtils;
93 titrationUtils = new TitrationUtils(protonatedAssembly.getForceField());
94 titrationUtils.setRotamerPhBias(298.15, pH);
95 for (Residue residue : protonatedAssembly.getResidueList()) {
96 String resName = residue.getName();
97 if (resNumberList.contains(residue.getResidueNumber())) {
98 if (resName.equalsIgnoreCase("ASH") || resName.equalsIgnoreCase("GLH")
99 || resName.equalsIgnoreCase("LYS") || resName.equalsIgnoreCase("HIS")
100 || resName.equalsIgnoreCase("CYS")) {
101 residue.setTitrationUtils(titrationUtils);
102 }
103 }
104 }
105 if (potentialEnergy instanceof OpenMMEnergy openMMEnergy) {
106 boolean updateBondedTerms = forceField.getBoolean("TITRATION_UPDATE_BONDED_TERMS", true);
107 openMMEnergy.getSystem().setUpdateBondedTerms(updateBondedTerms);
108 }
109 potentialEnergy.energy();
110 return protonatedAssembly;
111 }
112
113 public MolecularAssembly[] getProtonatedAssemblies() {
114 logger.info("Getting protonated assemblies");
115 MolecularAssembly molecularAssembly = getProtonatedAssembly();
116 List<Character> altLocs = protonFilter.getAltLocs();
117 int locs = 1;
118 if(altLocs!=null){
119 locs = altLocs.size();
120 for (int i = 0; i < locs; i++) {
121 if (altLocs.get(i) == null) {
122 altLocs.remove(altLocs.get(i));
123 }
124 }
125 }
126 MolecularAssembly[] molecularAssemblies = new MolecularAssembly[locs];
127 molecularAssemblies[0] = molecularAssembly;
128 for (int i = 0; i < altLocs.size(); i++) {
129 if (i != 0) {
130 logger.info(filename);
131 MolecularAssembly newAssembly = new MolecularAssembly(filename);
132 newAssembly.setForceField(forceField);
133 File structureFile = new File(filename);
134 protonFilter = new PDBFilter(structureFile, newAssembly, forceField,
135 forceField.getProperties(), resNumberList);
136 logger.info(newAssembly.getResidueList().toString());
137 protonFilter.setRotamerTitration(true);
138 protonFilter.setAltID(newAssembly, altLocs.get(i));
139 protonFilter.readFile();
140 logger.info(newAssembly.getResidueList().get(0).getAtomList().get(0).getAltLoc().toString());
141 protonFilter.applyAtomProperties();
142 newAssembly.finalize(true, forceField);
143 potentialEnergy = ForceFieldEnergy.energyFactory(newAssembly);
144
145 TitrationUtils titrationUtils;
146 titrationUtils = new TitrationUtils(molecularAssembly.getForceField());
147 titrationUtils.setRotamerPhBias(298.15, pH);
148 for (Residue residue : molecularAssembly.getResidueList()) {
149 String resName = residue.getName();
150 if (resNumberList.contains(residue.getResidueNumber())) {
151 if (resName.equalsIgnoreCase("ASH") || resName.equalsIgnoreCase("GLH")
152 || resName.equalsIgnoreCase("LYS") || resName.equalsIgnoreCase("HIS")
153 || resName.equalsIgnoreCase("CYS")) {
154 residue.setTitrationUtils(titrationUtils);
155 }
156 }
157 }
158 potentialEnergy.energy();
159 molecularAssemblies[i] = newAssembly;
160 }
161 }
162 return molecularAssemblies;
163 }
164
165
166 public boolean excludeExcessAtoms(Set<Atom> excludeAtoms, int[] optimalRotamers,
167 List<Residue> residueList) {
168 boolean isTitrating = false;
169 int i = 0;
170 for (Residue residue : residueList) {
171 Rotamer rotamer = residue.getRotamers()[optimalRotamers[i++]];
172 applyRotamer(residue, rotamer);
173 if (rotamer.isTitrating) {
174 isTitrating = true;
175 AminoAcidUtils.AminoAcid3 aa3 = rotamer.aminoAcid3;
176 residue.setName(aa3.name());
177 switch (aa3) {
178 case HID, GLU -> {
179
180 Atom HE2 = residue.getAtomByName("HE2", true);
181 excludeAtoms.add(HE2);
182 }
183 case HIE -> {
184
185 Atom HD1 = residue.getAtomByName("HD1", true);
186 excludeAtoms.add(HD1);
187 }
188 case ASP -> {
189
190 Atom HD2 = residue.getAtomByName("HD2", true);
191 excludeAtoms.add(HD2);
192 }
193 case LYD -> {
194
195 Atom HZ3 = residue.getAtomByName("HZ3", true);
196 excludeAtoms.add(HZ3);
197 }
198 case CYD -> {
199
200 Atom HG = residue.getAtomByName("HG", true);
201 excludeAtoms.add(HG);
202 }
203 default -> {
204 }
205
206 }
207 }
208 }
209
210 return isTitrating;
211 }
212
213
214 }