Build> delete atoms ///WOf course, in real life, water molecules are important, and we should treat them with more respect. Next we automatically add any missing bonds, based on interatomic distances: if two atoms are close enough that the distance between them looks like a bonding distance, the program will form a bond between them, if they are not already bonded, and guess a bond order for the new bond:
Build> set autobond mode chainFrom the messages in the message box we can see that only one bond was added (the one that we knew was missing). The set autobond mode chain command was crucial: without it the program would have added many more intermolecular bonds between the two inhibitor molecules /4apr/A/I and /4apr/A/X (the autobond method is not chemically intelligent...). If the autobond mode is set to chain automatic bonds are added only between atoms in the same chain (i.e. if they are in residues which have the same "parent" in the molecular data tree). After adding automatic bonds, we set the autobond mode back to its default (object, i.e. add bonds between atoms which are subject to the same display transformation).
Build> bond automatic *
Build> set autobond mode object
Hydrogens can be added in a semi automatic fashion, which usually works well, provided that the bond orders have been assigned correctly:
Build> modifyThe addition of hydrogens is semi automatic, since we had to give the program hints about how to add H's to groups with ambiguous protonation states. The hints can be given either by explicitly defining the formal charges, or, as we have done in the current example, by defining the ligance (i.e. number of bonds) of the ambiguous atoms. E.g., we have defined the ligance of arginine sidechain nitrogen atoms arg_*/N%* to be 3. The current number of bonds for arginine NE is 2 (so 1 hydrogen will be added to reach the requested ligance of 3), the current number of bonds to arginine NH1 and NH2 is 1 (so 2 hydrogens will be added to each to reach the requested ligance of 3). The atom selections in the atom ligance command have made use of the % wildcard character, which matches any single character (so the mainchain oxygen atoms O of Asp and Glu O are not included in * -atom name OXT asp_*,glu_*/O%* -done, while they would have been included had we used the selection * -atom name OXT asp_*,glu_*/O* -done (the * wildcard matches any sequence (including sequences of zero length) of characters. Well, it would not have made a difference, since the mainchain O's happen to have ligance 1 too, but I thought it would be a nice opportunity to introduce you to the % wildcard). The command sequence above produces charged N- and C-termini. In real life we would now have a look at the histidine sidechains, and try to guess their probable protonation state from their environment in the structure. To simplify things for the purposes of this tutorial we leave histidines at their default state (neutral, H on NE).
what: atom ligance 1 * -atom name OXT asp_*,glu_*/O%* -done
what: atom ligance 3 * -atom name arg_*/N%* -done
what: atom ligance 4 * -atom name lys_*/NZ ala_1/N -done
what: hydrogens *
what: done
The hydrogens command added H's with their default name (H). This would be o.k. for the energy minimization that we are going to do in the next section. However, some programs are a bit peculiar about non-unique atom names, so it is a good idea to automatically assign more sensible names:
what: atom name * -done -automatic hheavHydrogens will be named after the heavy atom to which they are attached, adding an additional digit (1,2,3) as needed to make names of H attached to the same atom unique (e.g. the three H's on CG1 of valine will be renamed to HG11, HG12, HG13, the H's on CA will be renamed to HA (HA1, HA2 for glycines). This will produce names for hydrogens which conform to the old PDB conventions.
Next, we are going to assign force field atom types:
what: atom type automatic taff * -doneCheck the message box for complaints about atoms which could not be assign automatically. If you have followed the example without making mistakes, all atoms should be correctly assigned now.
type: done
And finally we assign partial atomic charges by
what: atom q * -done mpeoei.e. partial atomic charges (q) are computed using the MPEOE method, which is tuned to produce charges compatible with the CHARMm 22 force field. From the message box, we can learn that the whole system has a total charge of -10 (formal and partial charges are identical, as they should be, if the method used to estimated the charge is of any value). The total charge is rather large, and in serious modeling, this would be cause for some thought (e.g. protonation state of the two catalytic Asp's, or of the His sidechain of the inhibitor, placement of counter ions, etc...), but in the spirit of HTM (high throughput modeling) we are happy to ignore these complications.
Now, we remove atom labels from the display (which have been put there by the modify commands) and return to the main Wit!P menu:
what: label remove *Now, this was hard work, and it would be a shame if we would lose it. So, let's save our system.
what: done
Build> end
Wit!P> freeze /The write command will store real coordinates, and therefore we freeze / (all display transformations, "object" and "world") to turn current display coordinates into real coordinates, so that when we read in the structure again it will already be centered and we do not have to use display find again.
Wit!P> write mol2 4apr.mol /4apr
We are using the mol2 file format, which is able to correctly store bond types and atom types partial atomic charges (the essential information needed for the energy minimization tat we are going to do next). However, the mol2 format can handle only up to four levels of the molecular data tree, so we will lose the /4apr/A level. Which is o.k., since it serves a crystallographic purpose only, and we won't need the crystal information (space group, unit cell geometry,...) any more.
If you see the molecule /4apr written to 4apr.mol message, it is save to exit Wit!P for a while, and have a break:
Wit!P> exit yes
We will do energy minimization after the break. Stay tuned...