<HTML><BODY style="word-wrap: break-word; -khtml-nbsp-mode: space; -khtml-line-break: after-white-space; "><DIV><DIV>On 16.07.2007, at 18:25, Jaroslaw Kalinowski wrote:</DIV><BR class="Apple-interchange-newline"><BLOCKQUOTE type="cite"></BLOCKQUOTE></DIV><DIV><BLOCKQUOTE type="cite"><DIV>programs to  calculate Poisson-Boltzman electrostatics), second it is</DIV><DIV>sometimes easier to use plain arrays rather than objects. For example</DIV><DIV>ParticleProperties are fine for representing properties defined for each</DIV><DIV>atom in a universe, but how could one represent data for selected objects</DIV><DIV>in the universe, or how could one represent two-atom properties like bond</DIV><DIV>orders defined for each pair? Arrays seem to be a choice, but for that one</DIV></BLOCKQUOTE><DIV><BR class="khtml-block-placeholder"></DIV>ParticleProperties objects behave like a dictionary, so I'd say the first choice for subsets of the whole system would be a dictionary mapping atom objects to the associated property value. For atom pairs, I'd use a tuple of atoms as the key.</DIV><DIV><BR class="khtml-block-placeholder"></DIV><DIV>Of course, there is nothing wrong with using arrays and atom indices, as long as you keep in mind that every time your program is run, the atom indices could be assigned differently in general. Note also that atom indices are not defined before the first call to universe.configuration(); in case of doubt, call it explicitly.</DIV><DIV><BR><BLOCKQUOTE type="cite"></BLOCKQUOTE><BLOCKQUOTE type="cite"><DIV>1. Ordering given by u.atomList() is not the same as ordering of</DIV><DIV>configuration array, is it?</DIV></BLOCKQUOTE><DIV><BR class="khtml-block-placeholder"></DIV>No.</DIV><DIV><BR><BLOCKQUOTE type="cite"><DIV>2. If not than how one can iterate over atoms in the same order as is used</DIV><DIV>for arrays?</DIV></BLOCKQUOTE><DIV><BR class="khtml-block-placeholder"></DIV>sorted_atoms = copy(universe.atomList())</DIV><DIV>sorted_atoms.sort(lambda a, b: cmp(a.index, b.index))</DIV><DIV>for atom in sorted_atoms:</DIV><DIV>    print atom</DIV><DIV><BR><BLOCKQUOTE type="cite"><DIV>3. Is it possible to get atom number in pdb file for peptide chains after</DIV><DIV>they have been created using PDBPeptideChain.createPeptideChain?</DIV></BLOCKQUOTE><DIV><BR class="khtml-block-placeholder"></DIV>No; in fact, there is no guarantee that every atom even has a PDB serial number, e.g. if the MMTK model has hydrogens but the PDB file hasn't.</DIV><DIV><BR class="khtml-block-placeholder"></DIV><DIV>There is one particular situation in which it is possible to access the PDB atom numbers: if </DIV><DIV><BR class="khtml-block-placeholder"></DIV><DIV>1) you use PDBConfiguration.createAll()</DIV><DIV>2) the PDB file contains all hydrogen atoms of the MMTK model</DIV><DIV>3) you don't add anything else to the universe</DIV><DIV><BR class="khtml-block-placeholder"></DIV><DIV>then the internal atom indices are in the same order as the atoms in the PDB file. This feature is used for interfacing with other simulation programs, e.g. when reading DCD trajectory files.</DIV><DIV><BR><BLOCKQUOTE type="cite"><DIV>4. Does pickling and unpickling preserve atomList ordering and/or array</DIV><DIV>ordering?</DIV></BLOCKQUOTE><DIV><BR class="khtml-block-placeholder"></DIV>Both. Note however that nothing guarantees that atomList ordering remains constant even between two successive calls to atomList(). It does in the current implementation, but it might not do so in the future.</DIV><DIV><BR><BLOCKQUOTE type="cite"><DIV>5. Does using copy.deepcopy do?</DIV></BLOCKQUOTE><DIV><BR class="khtml-block-placeholder"></DIV></DIV><DIV>If you copy the whole universe, it does in practice although no written rule guarantees that it does. It might not do so in the future.</DIV><DIV><BR><BLOCKQUOTE type="cite"><DIV>6. Does saving using PDBOutputFile and than reading using PDBConfiguration do?</DIV></BLOCKQUOTE><DIV><BR class="khtml-block-placeholder"></DIV>No.</DIV><DIV><BR><BLOCKQUOTE type="cite"><DIV>7. And can one know in advance what would be atom ordering in the pdb</DIV><DIV>written?</DIV></BLOCKQUOTE><DIV><BR class="khtml-block-placeholder"></DIV>No.</DIV><DIV><BR><BLOCKQUOTE type="cite"><DIV>8. What is the atom ordering in NetCDF trajectories?</DIV></BLOCKQUOTE><DIV><BR class="khtml-block-placeholder"></DIV>It is given by the system description in the netCDF variable "description". Every atom object in that description contains the index of that atom in the netCDF arrays.</DIV><DIV><BR><BLOCKQUOTE type="cite"><DIV>9. Can one assume that each residue has a continuous slice in</DIV><DIV>configuration array?</DIV></BLOCKQUOTE><DIV><BR class="khtml-block-placeholder"></DIV>No, except in the special situation described above.</DIV><DIV><BR><BLOCKQUOTE type="cite"><DIV>10. If so, than is the atom ordering in such slice identical for each</DIV><DIV>residue of a given type?</DIV></BLOCKQUOTE><DIV><BR class="khtml-block-placeholder"></DIV>Only of the PDB file has identical atom order for all residues of a given type.</DIV><DIV><BR class="khtml-block-placeholder"></DIV><DIV><BR><BLOCKQUOTE type="cite"><DIV>Some more problems.</DIV><DIV><BR></DIV><DIV>Is it possible to store additional data in a trajectory? I have a QM</DIV><DIV>forcefield that produces atomic charges, bond orders etc. Is it possible</DIV><DIV>to have my function called to provide data for each trajectory frame?</DIV></BLOCKQUOTE><DIV><BR class="khtml-block-placeholder"></DIV>The integrators and minimizers in MMTK don't do this at the moment. It should be possible to add such a a feature by writing a suitable action object, but I'd better don't promise this without some investigation.</DIV><DIV><BR class="khtml-block-placeholder"></DIV><DIV>If you write trajectories step by step from Python (using the SnapshotGenerator()), you can add any data you wish.</DIV><DIV><BR class="khtml-block-placeholder"></DIV><DIV><BLOCKQUOTE type="cite"><DIV>I would like to try doing QM/MM calculations. It seems that it is not</DIV><DIV>possible to apply Amber forcefield to the part of a system, so I was</DIV></BLOCKQUOTE><DIV><BR class="khtml-block-placeholder"></DIV>You can calculate energies and forces for a subset of the system. However, there is currently no way to tell the integrators to use this feature. If all you need is a simple velocity verlet algorithm, you could just rewrite it in a few lines of Python (it's fast enough) and calculate energies for subsets.</DIV><DIV><BR class="khtml-block-placeholder"></DIV><DIV>The "proper" way to do QM/MM in MMTK is to write a new force field. That force field would exclude the QM subsystem from the Amber force field term list and add the QM energies and forces instead.</DIV><DIV><BR><BLOCKQUOTE type="cite"><DIV>thinking about cloning MM part and adding QM part to a copy. Than I could</DIV><DIV>just write a forcefield that would just copy energy and gradients from MM</DIV><DIV>universe and add my QM terms. Do you think it is a good approach?</DIV></BLOCKQUOTE><DIV><BR class="khtml-block-placeholder"></DIV>That sounds doable, but you have to make sure that this really gives you the QM/MM mixture that you want. The forces on the MM subsystem would include interactions with the QM subsystem, for example.</DIV><DIV><BR><BLOCKQUOTE type="cite"><DIV>It is not possible to change protonation state after PeptideChain has been</DIV><DIV>created, is it?</DIV></BLOCKQUOTE><DIV><BR class="khtml-block-placeholder"></DIV>No. As far as MMTK is concerned, different protonation states for a residue are different residues.</DIV><DIV><BR class="khtml-block-placeholder"></DIV><DIV>It is possible with moderate effort to replace a residue in a chain by another residue, and thus change protonation states. You could look at the routine that creates S-S bridges as an example. However, any operation that changes the number of atoms invalidates the atom indices. They will be reassigned when needed, but this takes time, as does the subsequent recalculation of the force field terms etc. Changing the system is thus an expensive operation.</DIV><DIV><BR class="khtml-block-placeholder"></DIV><DIV>If you need changing protonation states in the course of a simulation, it is easier and more efficient to start with the fully protonated state and modify the force field terms to make unwanted hydrogens disappear from all interactions.</DIV><DIV><BR><BLOCKQUOTE type="cite"></BLOCKQUOTE></DIV><DIV>Konrad.</DIV><DIV><SPAN class="Apple-style-span" style="border-collapse: separate; border-spacing: 0px 0px; color: rgb(0, 0, 0); font-family: Arial; font-size: 16px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; text-align: auto; -khtml-text-decorations-in-effect: none; text-indent: 0px; -apple-text-size-adjust: auto; text-transform: none; orphans: 2; white-space: normal; widows: 2; word-spacing: 0px; "><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">--</DIV><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">---------------------------------------------------------------------</DIV><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">Konrad Hinsen</DIV><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">Centre de Biophysique Moléculaire, CNRS Orléans</DIV><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">Synchrotron Soleil - Division Expériences</DIV><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">Saint Aubin - BP 48</DIV><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">91192 Gif sur Yvette Cedex, France</DIV><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">Tel. +33-1 69 35 97 15</DIV><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">E-Mail: <A href="mailto:hinsen@cnrs-orleans.fr">hinsen@cnrs-orleans.fr</A></DIV><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">---------------------------------------------------------------------</DIV><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; "><BR class="khtml-block-placeholder"></DIV><BR class="Apple-interchange-newline"></SPAN> </DIV><BR></BODY></HTML>