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