[交流] 【原创】MAPS中构筑lammps的无定形结构的python程序已有3人参与

#!/usr/bin/python
###########################################################################
#                        amorphousLAMMPS.py  -  description
#                        ----------------------------------
#   begin                : Tue Oct  7 18:23:23 EEST 2008
#   copyright            : (C) 2008 by Scienomics
#   email                : loukas.peristeras@scienomics.com
###########################################################################
#
###########################################################################
#                                                                         #
#   This program and all subroutines,  data, and  files used by it are    #
#   protected by copyright and hence may not be used, copied, modified    #
#   transmitted, inspected, or executed by any means including the use    #
#   of electronic data processing equipment,  xerography, or any other    #
#   methods  without the express  written permission  of the copyright    #
#   holder.                                                               #
#                                                                         #
#   Copyright (C) 2008 Scienomics S. A.                                   #
#                                                                         #
###########################################################################
#
#   $Id: amorphousLAMMPS.py,v 1.4 2010/01/30 08:32:20 lperi Exp $
#   $Log: amorphousLAMMPS.py,v $
#   Revision 1.4  2010/01/30 08:32:20  lperi
#   Fix version handling.
#
#   Revision 1.3  2010/01/26 22:12:53  lperi
#   Fix bug 3182.
#
#   Revision 1.2  2010/01/04 19:29:03  lperi
#   Update with recoil growth.
#
#   Revision 1.1  2008/10/27 08:37:28  jrh
#   Created
#
#   $Revision: 1.4 $
#   $Release$

# Script to demonstrate the Amorphous Builder functionality in
# combination with LAMMPS optimizer.
# - builds a 100mer chain (Ethylene monomer is used with Dreiding force
#   field assigned)
# - builds a amorphous system (2 chains, 0.79 g/cm3) is build
# - perform LAMMPS optimization on the resulting system

import os
import time
import Maps
import MapsChemistryDataModel
import MapsViewer
import MapsLAMMPS
import MapsAmorphousBuilder
from PyQt4 import *

class AmorphousBuilderProtocol():
  """ The protocol used to prepare an amorphous cell of a given
      set of molecules.
  """
  plugin=MapsAmorphousBuilder.AmorphousBuilderPlugin.get()
  builder=None

  def create(self, molecule, nchain, density, temperature):
    self.builder=MapsAmorphousBuilder.AmorphousBuilder(app)
    self.builder.setFlexible(True)
    self.builder.setSimultaneously(True)
    self.builder.setNewProject(False)
    self.builder.setTemperature(temperature)
    self.builder.setDensity(density)
    self.builder.setInsertLateral(False)
    self.builder.setNumberOfInsertionTrials(50)
    self.builder.setNumberOfBuildTrials(5)
    self.builder.setNumberOfMoleculeBuildTrials(5)
    self.builder.setNumberOfBackSteps(5)
    self.builder.setNumberOfTrials(12)
    self.builder.setRecoiGrowth(True)
    self.builder.setScaleGroup(1.0)
    self.builder.setScaleLateralBonds(1.0)
    self.builder.setScaleBonds(1.0)
    self.builder.useNonBondedInteractions(True)
    self.builder.initializeRandomNumber(1234);

    molecules = [molecule]
    nmolecules = [nchain]
    self.builder.setMolecules(molecules,nmolecules)

    self.plugin.setBuilder(self.builder)
    self.builder.runningInThread=False

    app.statusBar().showMessage("Building amorphous system...")
    QtGui.qApp.processEvents()
    self.builder.startBuilder(False)
    QtGui.qApp.processEvents()
    app.statusBar().showMessage("Ready")
    QtGui.qApp.processEvents()
    if ( self.builder.getMolecule() ):
      self.builder.getMolecule().setName("amorphous")
      return(1)
    else:
      return(0)

  def getMolecule(self):
    if self.builder:
      return(self.builder.getMolecule())
    else:
      return(None)

class OptimizationProtocol:
  """ The protocol used to optimize the system.
  """
  plugin=MapsLAMMPS.LAMMPSPlugin.get()
  results=None
  molecule=None
  jobName=""
  projectName=""

  def setup(self):
    writer=MapsLAMMPS.LAMMPSWriter(app, self.plugin)
    writer.setOption("coordinateFile", self.jobName + ".data")
    writer.setOption("cutoff", "11.0")
    writer.setOption("dielectric", "1.0")
    writer.setOption("energyConvergence", "0.0001")
    writer.setOption("ensemble", "NVE")
    writer.setOption("forceCalculations", "5000")
    writer.setOption("forceConvergence", "0.000001")
    writer.setOption("forceFile", self.jobName + ".frcdump")
    writer.setOption("initialTemperature", "298.15")
    writer.setOption("method", "conjugatedGradient")
    writer.setOption("restartAtEnd", self.jobName)
    writer.setOption("scale12", "0.0")
    writer.setOption("scale12Coulomb", "0.0")
    writer.setOption("scale13", "0.0")
    writer.setOption("scale13Coulomb", "0.0")
    writer.setOption("scale14", "1.0")
    writer.setOption("scale14Coulomb", "1.0")
    writer.setOption("seed", str(int(time.time())))
    writer.setOption("steps", "500")
    writer.setOption("switchDistance", "8.0")
    writer.setOption("switching", "on")
    writer.setOption("timestep", "1")
    writer.setOption("title", self.jobName)
    writer.setOption("trajectoryFile", self.jobName + ".dump")
    writer.setOption("trajectoryFrequency", "100")
    writer.setOption("type", "forcefield")
    writer.setOption("typeOfCalculation", "optimization")
    writer.setOption("velocityFile", self.jobName + ".veldump")
    success=writer.createInput(self.results,self.molecule, self.jobName)
    if not success:
      QMessageBox.critical(app, "Note",  writer.getErrorDescription())
    return(success)

  def checkStatus(self):
    file=open("status", "r")
    line=file.readline()
    file.close()
    if line == "Completed":
      return(1)
    else:
      return(0)

  def optimize(self,molecule):
    self.molecule=molecule
    self.jobName=molecule.getName()+"_opt"
    self.projectName=cdm.findProjectNameForMolecule(self.molecule)
    # Add an experiment
    project=cdm.findProjectByName(self.projectName)
    experiment=Maps.Experiment(cdm,None,project,DataObject.END_TYPE)
    experiment.setName("LAMMPS " + str(self.molecule.getName()) + " optimization")
    project.addExperiment(experiment)
    self.results=Maps.ExperimentalResult(cdm,experiment)
    self.results.setName("LAMMPS results for " + str(self.molecule.getName()))
    experiment.addExperimentalResult(self.results)
    app.statusBar().showMessage("Optimizing structure using LAMMPS ...")
    # Setup LAMMPS optimization
    dirName=str(workingDirectory + "/" + self.jobName)
    if not os.path.exists(dirName):
      os.mkdir(dirName)
    os.chdir(dirName)
    success=0
    if self.setup():
      if os.name == "posix":
        prefix="./"
      else:
        prefix=""
      os.system(prefix + "run.py")
      if self.checkStatus():
        m=self.plugin.loadFile(self.results, \
                             self.jobName.toLatin1() + ".log")
        success=1
      else:
        QtGui.QMessageBox.critical(app, "Note", "The optimization using LAMMPS failed.")
    os.chdir("..")
    app.statusBar().showMessage("Ready")
    QtGui.qApp.processEvents()
    return(success)

def main():
  # Create output project
  project=MapsChemistryDataModel.Project(cdm)
  project.setName("Python polymer builder")
  cdm.addProject(project)

  # Find path to MAPS installation
  for d in app.getInstallPath():
    dir=d.toLatin1()
    if os.path.exists(dir): break
  file = "Ethylene.cml"
  fileName=dir + "/groups/Monomers/" + file

  # Build 100 mer chain
  print "Building polymer ..."
  molecule=MapsChemistryDataModel.Molecule(cdm)
  project.addMolecule(molecule)
  viewerPlugin.start(molecule)
  viewer=viewerPlugin.getActiveViewer()
  QtGui.qApp.processEvents()
  chainLength=100
  molecule.setName(str(chainLength) + "mer of " + file[0:len(file)-4])
  monomer=MapsViewer.Monomer(fileName, "Monomer", chainLength)
  monomers=[]
  monomers.append(monomer)
  viewer.buildPolymer(monomers, 1, 1, 0, MapsViewer.Viewer.ISOTACTIC)
  del monomers
  QtGui.qApp.processEvents()

  # Assign Dreiding force field
  print "Assigning force field ..."
  ff=MapsChemistryDataModel.ForceField(cdm, dir + "/forcefields/Dreiding.params")
  handler=MapsChemistryDataModel.ForceFieldHandler(app,ff.getForceFieldPath(),"")
  nullarg=None
  handler.setMolecule(molecule,nullarg)
  atomType=True
  atomCharge=True
  if handler.assignForceField(atomType,atomCharge):
    QtGui.QMessageBox.information(app, "Note", "The Dreiding force field was assigned.")
  else:
    QtGui.QMessageBox.warning(app, "Note", "There was a problem with the assignment\nof the force field.")
  QtGui.qApp.processEvents()

  # Create a amorphous system of 5 chains (0.79 g/cm3)
  nchain=2
  density=0.79
  temperature=450.0

  builder=AmorphousBuilderProtocol()
  print "Building amorphous system ..."
  if not builder.create(molecule,nchain,density,temperature):
    QtGui.QMessageBox.information(app, "Note", "Amorphous builder failed.")
  else:
    print "Optimizing amorphous system ..."
    optimizer=OptimizationProtocol()
    if not optimizer.optimize(builder.getMolecule()):
      QtGui.QMessageBox.information(app, "Note", "Optimization failed.")
  QtGui.qApp.processEvents()

#-------------------------
if __name__ == '__main__':
  # Check if MAPS is at least at version 3.2
  try:
    version=Maps.MapsApp.getVersion()
  except (TypeError, AttributeError):
    version=0
  if version < 302:
    QtGui.QMessageBox.information(None, "Note", "This script requires at least MAPS version 3.2.")
  else:
    # Get plugins needed
    app=Maps.MapsApp.getContainer()
    cdm=MapsChemistryDataModel.ChemistryDataModel.get()
    viewerPlugin=MapsViewer.ViewerPlugin.get()
    workingDirectory=app.getPreference("Working directory")
    main()

Originally posted by zyj8119 at 2010-12-26 17:41:04:
[code]#!/usr/bin/python
###########################################################################
#                        amorphousLAMMPS.py  -  description
#                        --------- ...

Originally posted by amynihao at 2010-12-28 09:11:21:

专家您好，这个程序是用来产生无定形结构的吗？那硅酸盐的玻璃态的结构可以用这个程序吗？这个程序是lammps里的，可以拟合到MS里吗？