From amasunov-0at0-shiva.Hunter.CUNY.EDU Tue Feb 20 22:48:28 1996 Received: from hcrelay.hunter.cuny.edu for amasunov <-at-> shiva.Hunter.CUNY.EDU by www.ccl.net (8.7.1/950822.1) id VAA04538; Tue, 20 Feb 1996 21:52:18 -0500 (EST) Received: from shiva.hunter.cuny.edu (amasunov(-(at)-)shiva.hunter.cuny.edu [146.95.128.96]) by hcrelay.hunter.cuny.edu (8.6.12/george0995) with SMTP id VAA10113; Tue, 20 Feb 1996 21:56:39 -0500 Date: Tue, 20 Feb 1996 21:53:07 -0500 (EST) From: Artem Masunov To: Computational Chemistry List cc: Artem Masunov Subject: Summary: CRYSTAL & all Message-ID: Dear Netters, Many thanks for your valuable replies. I tried to collect them all and added several URL sites, so this summary is rather long. I tried to edit, remove repititions, etc. but still... -------------------- Original Question -------------------- Recent summary on CRYSTAL & EMBED rises the question: are there any other programs available for periodical systems (I mean more advanced treatment then just empirical force field)? Unfortunately, Crystal92 does not optimize geometry, and Crystal95 is available only for mysterious EEC Human Capital and Mobility Network members (as Crystal Homepage says). AMPAC 2.2 had an option for 1D-periodicals, but it disappeared in later versions instead of expansion into 2,3D-cases. Does anybody know something else? ------------ My Comments and Relevant Information -------------- Besides Crystal95 and EMBED, Tourino group was distributing two DFT codes: PW and LAPW(=WIEN) Recently in Honolulu (Pacifichem'95) N.Tajima, H.Hirano, S.Tsuzuki & K.Tanabe presented a poster entitled "Totally ab initio prediction of the structures of CO2 molecular crystal". Authors used their own software. PROMET by A.Gavezzotti is not on WWW, contact address could be found on URL http://www.wkap.nl/natopco/NL50MEET.htm Info on MCI's POLYMORPH: http://www.msi.com/support/materials/applic/polymorphAnal/c2_PMAnal.html Ab initio method for geometry optimization of translationally symmetrical finite clusters: http://www.chem.joensuu.fi/people/juha_muilu/Research/ab_initio_method.html Info on ADF_BAND is available on URL http://iodine.chem.vu.nl/adf.html Info on CASTEP is available on URL http://www.tcm.phy.cam.ac.uk/castep/ Info on CAMP is available on URL http://hackberry.chem.niu.edu/ChemistrySoftware/MolecularDynamicsPrograms/CAMP another (all-electron) implementation of Car-Pimentallo scheme is PAW: http://www.zurich.ibm.com/~blo/ Seems that not only plane vawes but any quantum molecular dynamic program (such as Argus (http://hackberry.chem.niu.edu/ChemistrySoftware/SemiempiricalPrograms/argus) or QMDCP (MOTECC91) is capable to optimize crystal structure from the first principles... ---------------- New Questions ---------------- 1. Can MOPAC93/MOPAC7 optimize 3D-periodicals? 2. What semiempirical/ab-initio MD software is available? 3. Does anybody know e-mail of H.Hirano, S.Tsuzuki or K.Tanabe? ------------------- Replies ------------------- From: Xiaoming Zheng unsw.EDU.AU> Why not use the Car-Parrinello approach and their alternatives? You can look at Review of Modern Physics 64(1992)1045 and reference therein. The software implementing this approach calls CAMP-Atami by Japanese computational physics group (cost $25-$50). From: Andreas Ehlers there is a good program called ADF_BAND that can calculate periodic structures using the density funtional theory. From: Ole Swang The group of Baerends in Amsterdam maintains a program (ADF_BAND) for calculation of 1d, 2d, and 3d band structures. It uses DFT and a combination of Slater and numerical basis functions, and gives quite good results. Unfortunately, it is rather resource-demanding. Contact Bert de Velde (tevelde-0at0-chem.vu.nl) for further information. From: Andrey Khavryutchenko MOPAC 6 has 1D peridicity capability. It seems MOPAC 7 got 2D slabs and 3D solids, but I'm not sure. From: Fransisco Carlos Lavarda Also for 1D periodic systems you have the Mopac Package, versions 6, 7, and 93. From: Edoardo Apra' CRYSTAL92 is a Hartree-Fock program for the study of crystalline systems; to get informations about it, send an e-mail to CRYSTA92 at.at itocsivm.csi.it ------------------------------------------------------------ From: Georg Ostertag I know a group at the university of Ulm. They are developing a semiempirical method. As I know, they are now programming a facility to calculate cristals or polymers. I think, they are also interested in having contact with people using their method or program. The people who are working on the program are Rainer Mallon (Dipl. Physiker) and Manfred Doser (student) E-Mail: Manfred.Doser "at@at" physik.uni-ulm.de manfred "-at-" model2an.physik.uni-ulm.de Rainer.Mallon(-(at)-)physik.uni-ulm.de From: In BIOSYM technologies we have a good ab-initio software for PBC system. DSOLID: DFT soft based on the numerical basis set for PBC systems Plane_Wave is a first-principles (ab initio) quantum mechanics soft package that performs accurate calculations on solid-state systems. Plane_Wave calculates variational self-consistent solutions to the density functional theory (DFT) equations for insulating and semiconducting systems that are subject to 3D periodic boundary conditions. The method uses a plane-wave expansion for the wavefunctions and a pseudopotential representation of the ions. ESOCS is a first-principle (ab-initio) quantum mechanics packages that performs accurate theoretical calculations on a wide range of 3D solid-state systems including metals, semiconductors and surfaces. ESOCS stands for Electronic Structure of Close-packed Solids. It is based on spin density functional methods and the Atomic Sphere Approximation (ASA). The distinguishing feature of the code is a spherical wave expansion of the wave functions, which are centered at each atomic state of the solids.Because of this construction, the projection of quantities like densities of state onto atomic sites is particularly easily computed. Since ESOCS does not use the frozen core approximation, changes in the core spectra due to core relaxation can be calculated. From: Lasse Hemmingsen Concerning your question about first principle quantum chemical methods on periodic systems (crystal-like that is), I suggest that you contact Prof. Karlheinz Schwarz or Peter Blaha (Prof. Schwarz's email adress is kschwarz &$at$& email.tuwien.ac.at). Their program, WIEN9?, is made for calculations on periodic systems, and it has given excellent results for example in calculating electric field gradients. From: heinz()at()olymp.theochem.tuwien.ac.at Reply to: Karlheinz Schwarz Our WIEN95 does allow a geometry optimization, since force are implemented into the program. The optimization, however, is presently tested and thus is not fully automated yet. Please note that we will organize a workshop for the use of WIEN95 (10-13 April 1996, in Vienna). For information on the full-potential linearized augmented plane wave (LAPW) program WIEN95 see Worl Wide Web: http://www.tuwien.ac.at/theochem/wien95/wien95.html From: J Perlstein If you are interested in predicting the crystal structure of a molecule knowing only the atom-atom connectivity, there are at least four programs presently available for this. One called POLYMORPH from BIOSYM/MSI based on Karfunkels work another called PACK from Chemical Design, Inc. based on my own work. Both programs use a Monte Carlo simulated annealing algorithm to find the global and nearby local minima for a molecule in a particular space group. The former program should be adequate for rigid molecules, the latter allows for exocyclic conformational flexibility for molecules containing up to 12 torsion bonds. There is also a program available from D. Williams at the University of Kentucky available from QCPE (and I believe also as a commercial package from Williams) which uses energy minimization techniques and allows for some flexibility as well as some software from Angelo Gavezzotti at the University of Milan which uses a cluster approach to packing. PACK can handle molecules with up to 130 atoms and 12 exocyclic torsion bonds. It uses either MM2 for non-hydrogen bonded systems or CFF91 for hydrogen bonded systems with one molecule in the packing unit. PACK runs as an interface to the CHEMX/CHEMLIB molecular modeling program on either SGI's or IBM RS/6000. The CHEMX/CHEMLIB program is available from Chemical Design, Inc. in Mahwah,NJ. Contact Debbie Dunn at Chemical Design for further details.(Tel: (201)529-2443). There are several papers which discuss methods for predicting packing geometries of molecules in solids: 1)Gavezzotti, A. JACS 113, 4622 (1991) 2)JR Holden, Z. Du, HL Ammon, J. Comput. Chem. 14, 422 (1993) 3)HR Karfunkel, B. Rohde, FJJ Leusen, RJ Gdanitz, GJ Rihs, J.Comput.Chem 14, 1125 (1993) 4) J Perlstein, JACS 116, 11420, 1994 PACK is being used by a number of groups including Jenekhe for predicting packing of polymer chains (See Roberts et.al. Chem. Mater. 6, 658 (1994)) and the Whitten group who are interested in monolayer packing (See Chen et. al. Phys. Chem. 98, 5138 (1994)) From: "D. E. Williams" Computer program mpa (molecular packing analysis) predicts the structures of molecular clusters or crystals by energy minimization, based on a variety of published force fields or a user-provided force field. Molecular docking energy minimization and docking structure prediction is simply treated as a molecular cluster calculation where molecules in the cluster are not identical. The intermolecular or nonbonded energy of the molecular assembly is represented by a pairwise sum over atoms in different molecules. The program accepts either (exp-6-1) or (n-6-1) nonbonded inter- atomic potentials, referred to as Buckingham or Lennard-Jones functions. A torsional potential is accepted for rotations about selected internal bonds. The program can handle ionic species and provision is made for net atomic charges or lone pair electron site charges. The structural variables considered by the program are the rotations and translations of several molecules, and selected internal rotations. Molecules may be related by space group symmetry operations. For a crystal calculation lattice symmetry is added and the six (or fewer) parameters of the unit cell may be selected as variables. Crystal lattice energy calculations can be made in one, two, or three dimensions with any space group symmetry. There can be more than one independent molecule in the crystallographic asymmetric unit. Crystal lattice sums are accurately evaluated using the accelerated convergence method. Provision is made for summation over electrically neutral subunits. The crystal unit cell is checked for correctness before each contact table is generated. This involves tranformation of the structure to the conventional primitive reduced cell and moving the molecules as necessary to be inside the defined cell. First and second derivatives of the molecular cluster or crystal lattice energy are evaluated analytically (rather than numerically). The program uses several energy minimization methods, designed to find both local and global minima. Eigenvalues and eigenvectors of the second derivative matrix (hessian) are calculated. If all eigenvalues are positive (positive definite hessian) the Newton-Raphson method (NR) is used to find shifts toward the energy minimum. If the hessian is not positive definite off-ridge eigenvector minimization (OREM) is used. In this method energy is minimized along eigenvector directions corresponding to the most negative eigenvalues. In a few special situations steepest descents (SD) minimization is used. The program includes provisions for an annealing schedule. Annealing may be useful to carry a structure from a subsidiary energy minimum to its global minimum. In the annealing process the structure is given random shifts or thermal kicks which simulate the effects of thermal motion. The resulting higher energy structure is then relaxed, but before the energy is completely minimized additional random shifts are applied as specified by the annealing schedule. These thermal kicks occur at specified effective temperatures corresponding to a threshold energy decrement per minimization cycle. When the effective temperature decreases to zero the structure energy is fully minimized. Program mpg (molecular packing graphics) reads an output file from mpa to graphically display predicted molecular cluster or crystal structures. Mpg differs from the many currently available molecular graphics display programs in that it is specifically designed to show salient aspects of the packing structures of clusters and crystals in a convenient and informative manner. Limits on number of atoms, etc. The default maximums are as follows. These limits can be changed during program installation. atoms in input list 1000 molecules in asymmetric unit 64 space group operators 24 molecules in contact table 900 contacts in contact table 100000 neutral groups 200 number of variables 192 potential sets 30 internal rotations 10 Selected publications of Donald E. Williams D. E. Williams, "Molecular Packing Analysis", Acta Crystallographica 1972, A28, 629-635. D. E. Williams, "Optimally Spaced Rotational Grid Points", Acta Crystallographica 1973, A29, 408-414. D. E. Williams and T. L. Starr, "Calculation of the Crystal Structures of Hydrocarbons by Molecular Packing Analysis", Computers & Chemistry 1977, 1, 173-177. D. E. Williams and D. J. Houpt, "Fluorine Nonbonded Potential Parameters Derived from Crystalline Perfluorocarbons", Acta Crystallographica 1986, 286-295. D. E. Williams, "Alanyl Dipeptide Potential-Derived Net Atomic Charges and Bond Dipoles and their Variation with Molecular Conformation", Biopolymers 1990, 29, 1367-1386. D. E. Williams, "Science Citation Classic: Nonbonded Potential Parameters Derived from Crystalline Hydrocarbons", Current Contents 1990, 30(11) 14. D. E. Williams, "Net Atomic Charge and Multipole Models for the Ab Initio Molecular Electric Potential", Reviews in Computational Chemistry 1991, 2, 219-312. D. E. Williams, "OREMWA Prediction of the Structure of Benzene Clusters: Transition from Subsidiary to Global Energy Minima", Chemical Physics Letters 1992, 192, 538-543. D. E. Williams and T. R. Stouch, "Characterization of Force Fields for Lipid Molecules: Applications to Crystal Structures", Journal of Computational Chemistry 1993, 14, 1066-1076. D. E. Williams, "Accelerated Convergence of R(-n) Lattice Sums", International Tables for Crystallography 1993, Volume B, 374-382. Y. L. Xiao and D. E. Williams, "Genetic Algorithms for Docking of Actinomycin D and Deoxyguanosine Molecules with Comparison to the Crystal Structure of Actinomycin D-Deoxyguanosine Complex", Journal of Physical Chemistry 1994, 98, 7191-7200. D. E. Williams, "Ab Initio Molecular Packing Analysis", Acta Crystal- lographica 1996, Section A, in press. QUOTATION (for academic use only) mpa/mpg single computer license and installation . . . . . .$1,995.00 Installation requires telnet access (with temporary password) to the computer directory where the software is to be installed Test example calculations. Several test files for small molecules are available to illustrate various features of mpa. The first group of tests is for crystals starting from the observed structure. The second group illustrates crystal structure prediction, not starting from the observed structure. The third group is for prediction of molecular clusters and molecular docking. For each test there is a named input file (.mpa file extension) and three named output files (.out, .summary, and .restart). The .out file is the main output from the program, and the .summary and .restart files are renamed versions of mpa.summary and mpa.restart produced by the test. The computer time required for the examples will vary widely depending on the speed and capacity of the computer. The folowing times were obtained using a 100 MHz SGI Indy computer. Test Description Time(seconds) 1 carbon_dioxide 4 2 benzene 14 3 benzene_neutral 26 4 benzene_reciprocal 1 5 pyrimidine 48 6 hexafluorobenzene 175 7 triphenylbenzene 176 8 benzene_large_shifts 38 9 benzene_abinit 2670 10 urea_abinit 1595 11 benzene cluster 198 12 benzene global minimum 221 13 actinomycin docking 318 Index to the test examples. A. Crystal structure predictions starting from the observed structure. These calculations relax the structure to an energy minimum; they also test the accuracy of the force field selected. 1. Carbon dioxide 2. Benzene 3. Benzene: summation over neutral units 4. Benzene: reciprocal lattice sum 5. Pyrimidine: lone pair electron sites 6. Hexafluorobenzene: more than one molecule in the asymmetric unit 7. 1,3,5-triphenylbenzene: internal rotation B. Crystal structure predictions not starting from the observed structure. These calculations test the predictability of the crystal structure. Subsidiary energy minima, as well as a global minima, may be obtained. 8. Benzene, with large shifts 9. Benzene, ab initio prediction 10. Urea, ab initio prediction C. Molecular cluster predictions. 11. Benzene four molecule cluster 12. Benzene cluster global minimum 13. ActinomycinD-deoxyguanosine molecular docking 1. Carbon dioxide (carbon_dioxide.mpa, carbon_dioxide.out, carbon_dioxide.summary, and carbon_dioxide.restart) This is a quick and easy crystal structure calculation to see if the program is operational. The program gives a date and time stamp to identify the run. After the title the unmodified input atomic coordinates are listed followed by copies of the initial rotation matrix, initial rotation center, and initial translation. The center of mass is the new center of rotation but is not different from the initial in this case. The control parameters are then listed. There are three atoms in this molecule (NA=3) and four molecules in the unit cell (NS=4). Atomic coordinates are cartesian (NCEL=0) and no rotation or translation is specified (NROT=0). Net atomic charges are specified (NCHRG=1) and 40 contact tables are to be generated (NTAB=40). Carbon-hydrogen bonds are not to be foreshortened (XSHORT=0.0) and the initial energy is not known (EZERO=0.0). No annealing is specified (NAGITATE=0) and full output is requested (NPRINT=0). The specified force field potential energy parameters are then listed. Note that the user can supplement or customize force fields by editing the file mpa.pot. The program determines the number of molecules (rigid bodies, NTHE=1), the number of neutral subgroups (NGRP=1), and the number of subrotatations (rotations about bonds, NPSI=0). A check is made for net molecular electric charge. Carbon dioxide crystallizes in space group Pa3_, which has 24 symmetry operations. The reference molecule is placed at the origin pointing along the diagonal threefold axis of the cell; the molecular symmetry includes the crystallographic threefold and an inversion. When these 6 operations are deleted from the general operations (as shown in the International Tables for Crystallography), only three screw axis operations remain, plus the identity operation; the subgroup is equivalent to P212121 with a=b=c. The program notes that all symmetry matrices are diagonal (NSDIAG=0). When this is the case atom pair index IA will always be taken equal or greater than IB, and the energy and derivative values doubled when IA is not equal to IB. If any symmetry matrix has nonzero off-diagonal elements, NSDIAG is set to one, and all values of IA and IB are considered. The input lattice constants are printed, along with the lattice symmetry type and the cell volume. The center of mass of the molecule is found and output as the new rotation center (unchanged in this case). The first calculations conditions line requests energy minimization to end when the cycle decrement is less than 0.00001 (ZEND); contacts are to be considered up to 14A (SUMLIM). The accelerated convergence constants are 0.125 for coulombic interaction (CONV1) and 0.15 for dispersion interaction (CONV6). The reciprocal sum limit is 0.0 (RLIM) and the charge scale factor is 1.0 (SCALQ). The second calculations conditions line requests the use of cartesian coordinates (ICAR=0, rather than principal axes of inertia), summation over electrically neutral units (IREQ=1), no output of reciprocal lattice terms (IROUT=0), and final output of contacts less than 0.5A plus the sum of van der Waals radii (ILIST=1, ADJUST=0.5). The full hessian is to be calculated and output (NBLOCK=0) and this is not a molecular cluster calculation (NCLUSTER=0). If the hessian is not positive definite, one eigenvector corresponding to the most negative eigenvalue is to be used in OREM refine- ment (NEVNUM=1). Variable shifts are limited to 0.2A or rad. The "a" lattice constant (only) is selected. No rotations or translations are selected, since these are fixed by symmetry. Cell constants "b" and "c" are treated as dependent parameters, set equal to "a". by IDEP. The main calculation begins by scanning cells from -4 to +4 from the reference cell to generate a molecule list (LTSM=362 molecules) and a contact table (table 1, ITABW=2172 contacts). The initial lattice energy is -27.4392 kJ/mol (EZERO), which is broken down into a coulombic component of -15.9697 (ECOUL), a dispersion component of -39.4472 (EDISP), and a repulsion component of 27.9777 (EPACKR). The calculated values of the first (gradient) and second (hessian) derivatives are listed. In this simple case the eigenvalue of the hessian is just the second derivative itself. Since all eigenvalues are positive, a Newton-Raphson energy minimization cycle is initiated. The energy falls to -28.14341 by shifting the "a" lattice constant by 0.2000; note that the cubic symmetry is retained as specified by the dependent parameters, and that the shift was limited by XLIMIT. In cycle 2 the same contact table is used (to insure convergence, it is important that final refinement be based on same contact list) and only the first derivatives are calculated. Recalculation of the hessian is omitted because it is not expected to change very fast and its calculation is very time-consuming. The lattice constant decreases by 0.0625, showing the nonlinearity of the minimization and the usefulness of shift limits. In cycle 3 there is a further small adjustment of 0.0039 in the lattice constant. It is now time to generate a new contact table (ITAB=2) with 320 molecules and 1920 contacts, after the new center of mass is computed (in this case it is unchanged). Note that the initial energy (-28.2466) is slightly different from the ending energy of the previous cycle (-28.2464). The difference can be either positive or negative, and is a normal consequence of a change in the contact table. Smaller summation limits and larger variable shifts will increase the difference. In cycle 2 the variable shift becomes smaller than ZEND and the calculation terminates. The program lists the final energy values and the energy per molecule, and the molecular rotation and translation (here zero). A summary listing of short contacts is produced as requested by ILIST and ADJUST. Contacts are sorted in order of potential type. A summary of the lattice constant shifts, cell volume, and reduced cell is given. Of course, the lattice constant shifts reflect the accuracy of the force field. A list of the final atomic coordinates is given, both in cartesian and fractions of the unit cell. 2. Benzene (benzene.mpa) There are four molecules per cell in space group Pbca, with the molecules on inversion centers. The published structure lists only half of the atoms; the rest are generated by applying inversion symmetry. An idealized symmetrical planar molecule with averaged distances was fitted to the observed structure, and deletion of the inversion operation from Pbca yields the operators of P212121. Notice that the C-H distance has been set to a foreshortened 1.027A to be consistent with the WH86 force field. Alternatively, the C-H distance in the input coordinates could have been the normal 1.097A and XSHORT set to -0.070. Only rotations are allowed, since the molecule must be constrained to the inversion center. The angles of the orthorhombic cell are also not varied. The input molecule is rotated by a previously determined matrix which gives the best least squares fit to the observed orientation. The program considers 79 molecules with 2767 nonbonded contacts resulting in an initial energy of -51.6259. Since the eigenvalues of the hessian are all positive Newton-Raphson refinement is selected. In cycle 2 of the second contact table the energy drops less than 0.0002 and the energy minimization is finished. The final energy is -52.5376. The angular shift was 2.909 deg degrees from the input orientation. A list of short nonbonded contacts in the optimized structure is given. The a, b, and c lattice constants changed by 0.0676, -0.1134, and 0.3030A. The final cell and molecular volumes are 481.96 and 120.49 cubic A. The optimized atomic coordinates are listed in both cartesian and fractional cell coordinates. 3. Benzene, summation over neutral units (benzene_neutral.mpa). Using IREQ=1 and the same summation limit (10A), the program selects 80 (instead of 79) molecules which overlap a 10A sphere drawn around any reference molecule atom. There are now 6240 nonbonded contacts, requiring more time for the calculation. The summation over neutral units results in a more accurate initial energy of -51.5461. The final energy is -52.4576, the angular shift is 2.965 degrees, and the lattice constant shifts are 0.0697, -0.1139, and 0.3065. Setting NPRINT=1 in the command line reduces the amount of output. Setting NBLOCK=1 suppresses the printing of the hessian. 4. Benzene, reciprocal lattice sum (benzene_reciprocal.mpa). The reciprocal lattice sum is evaluated with the data of test 2. NTAB=0 causes the program to calculate the lattice energy only. The individual terms in the coulombic and dispersion reciprocal sums are listed (76 terms). The total reciprocal lattice sums are 0.0047 and -0.0036 to for coulombic and dispersion, respectively. This type of calculation is useful for selection of optimum values of the convergence constants conv1 and conv6. Generally, these constants should be chosen so as to make reciprocal lattice sum contributions neglegible relative to truncation error in the direct lattice sum. Note that mpa does not include reciprocal sum contributions to first or second derivatives; therefore conv1 and conv6 should be kept small. 5. Pyrimidine, with lone pair electron sites (pyrimidine.mpa) The molecule is in a general position in an orthorhombic lattice. Since the space group is polar along z this molecular translation is not selected. The atom list includes lone pair electron sites X1 and X3 0.25 A from the nitrogens. Only coulombic interaction is considered between the lone pair sites. Because the charges are large, IREQ is set to 1 to cause the program to sum over whole asymmetric units. This improves the accuracy of the coulombic sum. The calculation quickly converges with a molecular rotation of 1.719 deg and a molecular translation of 0.064 A. The a, b, and c cell constants change by 0.0007, 0.2556, and -0.0471 A. The energy at the input observed crystal structure is -56.5453 and the energy at the calculated structure has been minimized to -57.3145. 6. Hexafluorobenzene: more than one molecule in the asymmetric unit (hexafluorobenzene.mpa) Hexafluorobenzene crystallizes in P21/n with six molecules per cell. Since the space group has fourfold equivalent positions, one molecule is in a general position and the second molecule is on an inversion center. The mpa space group is P21 with three molecules in the input atom list (asymmetric unit). Molecule 1 and molecule 2 are related by inversion around one center, and molecule 3 sits on a second inversion center. To retain the observed symmetry, IDEP numbers are assigned. ISEL selects the monoclinic cell constants, the rotations and translations of molecule 1, and the rotations (only) of molecule 3. Molecule 2 is required by the IDEP numbers to rotate the same way as molecule 1, and to translate in the opposite way. This dependency retains the inversion relationship between molecules 1 and 2. By not allowing translation, molecule 3 is fixed on the other inversion center. Note that in the published x-ray structure only half of the atomic coordinates of molecule 3 are given; the others are easily generated. Similarly, all of the coordinates of molecule 2 are generated from the complete list of atomic coordinates for molecule 1. In the example file an idealized molecule was fitted to the observed structure to give the initial rotations and translations. 7. 1,3,5-Triphenylbenzene with intramolecular rotations (triphenbenz.mpa) This molecule is used to illustrate the use of subrotations, or internal rotations. The molecule is divided into four parts to enable summation over electrically neutral groups. Each of the three peripheral benzene rings is initially coplanar with the central ring. The nucleus central ring has atoms assigned IPSI=0 including three bonded carbons and three bonded hydrogens. Rotations are selected about the three phenyl-phenyl bonds. IPSI values are assigned 1, 2, or 3 for atoms in the peripheral rings. The specified subrotation potential has twofold periodicity with a maximum conjugation energy of -40.00 kJ/mol for coplanar phenyls. The observed cell and space group are used. Note that in this case EZERO values include intramolecular strain energy. Initially, the repulsion intramolecular energy ESTRNR is 208.2110, leading to EZERO=-24.7450. In table 4, ESTRNR decreases to 90.4128 and EZERO has decreased to -171.2733. Intramolecular coulombic and dispersion energy is also included in the calculation (ESTRN1 and ESTRN6). Including the shifts found from table 14, the final predicted phenyl rotation angles are 318.74, 38.33, and 319.06, in good agreement with the observed angles. Note one phenyl "propeller blade" is reversed in this structure. This structure can be used to illustrate the characteristics of summation over neutral groups to achieve greater accuracy in the coulombic energy summation. With four neutral groups specified, the contact table contains 38762 entries. If IREQ is set to 1 in the data file, then IGRP numbers are ignored and the contact table has 52854 entries. Compute time in the program is approximately proportional to the size of the contact table. If only an atom-atom SUMLIM is used (NGRP=1 and IREQ=0) there are only 8809 entries in the contact table. In this case the calculation is much faster, at the expense of poorer convergence of the coulombic sum because of summation over a non-neutral domain. 8. Benzene, with large shifts (benzene_large_shifts.mpa). Instead of the observed RMINP matrix, the program is given a unit matrix which gives a grossly incorrect molecular orientation with an initial energy of 123.9. The observed orthorhombic lattice constants and space group are initially specified. The cell constants and rotations of the molecule are varied. The initial hessian has four negative eigenvalues. OREM along eigenvector 1 reduces the energy to 27.96017, along eigenvector 2 to 10.73014, along eigenvector 3 to -9.53552, and along eigenvecto 4 to -12.89525. In contact table 2 the energy calculated from a new contact table is -13.1912. Three eigenvalues remain negative, and OREM decreases the energy to -38.66789. Contact table 3 results in two negative eigenvalues, and the energy decreases to -42.03437. Contact tables 4, 5 and 6 show single negative eigenvalues, and the energy decreases to -52.05996. A supplementary SD cycle is done because the energy didn't decrease much with the last table. With contact table 7 the calculation enters a NR region and quickly refines to a minimum of -52.5121. Of course, the rotation of the molecule is large (49.843 deg). The lattice constants agree with the relaxed observed structure. The conventional reduced primitive cell is given. In the .summary file the lines where Newton-Raphson refinement was used are appended with NR. 9. Benzene, ab initio molecular packing analysis space group prediction with Z=4 (benzene_abinit.mpa). In this test no observed information is used regarding unit cell constants, molecular position, or symmetry. To produce this test, first several runs were made with NROT=3, in which random initial angular orientation of all four molecules is specified. The initial cell is 10x10x10A and initially the molecules are in face centered positions. After several trials, the observed structure was predicted correctly. The Lattman angles from this successful trial were then inserted into the input file, with change of NROT to 5. There are 27 variables in this calculation, enumerated as follows: 6 cell constants, 3 rotations of the first molecule, and 6 rotations and translations of the subsequent molecules. The first molecule is not translated in order to set the origin of the lattice. The end of the .out file shows convergence to -209.5558, or -52.3890 per molecule. The cell angles are all 90 deg, and the cell edges agree with the observed values. To determine the space group, look at the final cell coordinate list. Noting that all C and H atoms are equivalent in benzene, molecule b is related to molecule a by -x, y+1/2, -z+1/2, which is a screw axis along y. Similar consider- ations lead to the identification of the observed molecular packing space group P212121. Now look at the .summary file. The initial energy for four molecules is -53.6005. At first, OREM shifts are made in six different directions. At table 5 the number of negative eigenvalues decreases to five. There are reduced cell changes at tables 7, 8, 11, 15, and 18. NR minimization kicks in at table 12, but is lost at table 14. NR kicks in again at table 15, only to be lost again at table 16. Beginning at table 17 a continuous NR path leads to the converged energy minimum. This calculation is further described in reference Williams (1996). 10. Urea, ab initio molecular packing analysis with Z=2 (urea_abinit.mpa). Since net atomic charges are large in this molecule, summation is over a neutral domain (IREQ=1). The Biosym force field is selected (IPOTR=4). As in the benzene example, a number of random trials were made with NROT=3 to locate minima. The lowest minimum found agrees with the observed structure with space group P4_21m. Note that this space group is not observed with high frequency in the Cambridge Structural Database of organic crystal structures. The test file uses NROT=5 and inserts successful trial Lattman angles. The starting model is body-centered in a 8x8x8A cell. There are 15 degrees of freedom and the initial energy is -16.6118 for two molecules. The summary file shows reduced cell changes at tables 4, 5, 9, and 14. The program switches to NR refinement at table 11, but reverts to OREM in table 13; at table 14 final NR refinement continues to an energy of -174.8567 for two molecules. The predicted tetragonal structure is in agreement with the observed structure; it can be displayed with mpg. This calculation is further described in reference Williams (1996). 11. Benzene 4-cluster subsidiary minimum (b4cluster.mpa) The starting model is tetrahedral, with the molecules on alternate corners of an 8A cube. The lattice constants are set to 200 A with the summation limit set to 20 A. A large shift limit of 1.0A is allowed. Rotation and translation of molecules 2, 3, and 4 are selected. There is only one lattice symmetry entry (LTSM=1) and 864 contacts (ITABW, (12x12)x6). The initial energy is barely negative, -0.5543. Us of up to 6 negative eigenvalues was specified, leading to EZERO=-1.5952 in table 2. At table 14 NR refinement begins, temporarily interrupted at tables 16 and 17. Final energy minimization at table 20 yields EZERO=-50.8925. This tetramer structure is subsidiary minimum (1b) of reference Williams (1992). 12. Benzene 4-cluster global minimum (b4anneal.mpa) The test is the same as test 11, but with annealing added (OREMWA) to assist in the search for the global energy minimum. The annealing schedule is started at 200K and reduced by 20K decrements, with initial maximum random shifts of 0.5A. Since annealing is a Monte Carlo process, the results of each run will be different. File b4anneal.out shows a run that successfully located the global minimum described in reference Williams(1992). Agitation begins immediately at 200K, leading to a starting energy of -0.5585. Energy minimization proceeds along the 6 selected eigenvectors, yielding -1.9849. Since the energy has not decreased by at least eagitate (5.9861), another agitation takes place at 180K. This process continues, lowering the annealing temperature, eventually leading to a T=0K model with an initial energy of -54.2866. This model quickly refines to the global energy minimum of -55.6281 reported in reference Williams(1992). 13. ActinomycinD-deoxyguanosine molecular docking (actinomycin.mpa) The starting model was taken from the crystal, but not including the waters of solvation. There are 176 atoms in the actinomycinD molecule and 32 atoms in each of the two deoxyguanosine molecules. The lattice constants are set to 200x200x200 so that adjacent cells are too far away to have any contacts included with SUMLIM=25. If SUMLIM is large enough, all atoms of each of the three molecules will be included in making up the contact table. In this case it can be verified that the contact table contains (176)(32)(2) + (32)(32) = 12288 entries. Since this is not a lattice, the convergence constants are set to zero. In addition to not varying the cell constants, actinomycinD rotations and translations are not varied. Each deoxyguanosine has six rotations and translations allowed, for a total of twelve variables. The starting energy is -55.2227. The calculation does not start in a NR region, as there are four negative eigenvalues. In table 2 the number of negative eigenvalues decreases to one, and in table 3 a NR region is reached. Refinement proceeds to a minimum energy of -129.4025. The two deoxyguanosine molecules have been rotated by 12.258 and 1.955 deg, and translated by 0.875 and 0.182 A. Of course, this model represents an isolated gas phase cluster and is not expected to agree exactly with the crystal configuration. In addition to deviations normally experienced because of approximations in the force field and neglect of thermal effects, waters of solvation are not included in this model. --------------------------- End of Replies -------------------------- Once again, thank you for replies. Artem __ _________ / \ / _ _ \ Artem Masunov - amasunov "-at-" shiva.hunter.cuny.edu / \\ \\ \\ \ Chemistry Department, Hunter College / /\ \\ \\ \\ \ City University of New York / ____ \\ \\ \\ \ 695 Park Avenue, New York, NY 10021 /__/\__/\__\\__\\__\\__\ Tel: (212) 725-0317, Fax: (212) 772-5332 \__\/ \/__//__//__//__/