39,59 €
Mehr erfahren.
- Herausgeber: Packt Publishing
- Kategorie: Wissenschaft und neue Technologien
- Sprache: Englisch
LAMMPS is one of the most widely used tools for running simulations for research in molecular dynamics. While the tool itself is fairly easy to use, more often than not you’ll need to customize it to meet your specific simulation requirements. Extending and Modifying LAMMPS bridges this learning gap and helps you achieve this by writing custom code to add new features to LAMMPS source code. Written by ardent supporters of LAMMPS, this practical guide will enable you to extend the capabilities of LAMMPS with the help of step-by-step explanations of essential concepts, practical examples, and self-assessment questions.
This LAMMPS book provides a hands-on approach to implementing associated methodologies that will get you up and running and productive in no time. You’ll begin with a short introduction to the internal mechanisms of LAMMPS, and gradually transition to an overview of the source code along with a tutorial on modifying it. As you advance, you’ll understand the structure, syntax, and organization of LAMMPS source code, and be able to write your own source code extensions to LAMMPS that implement features beyond the ones available in standard downloadable versions.
By the end of this book, you’ll have learned how to add your own extensions and modifications to the LAMMPS source code that can implement features that suit your simulation requirements.
Das E-Book können Sie in Legimi-Apps oder einer beliebigen App lesen, die das folgende Format unterstützen:
Seitenzahl: 378
Veröffentlichungsjahr: 2021
Ähnliche
Extending and Modifying LAMMPS
Writing Your Own Source Code
A pragmatic guide to extending LAMMPS as per custom simulation requirements
Dr. Shafat Mubin
Jichen Li
BIRMINGHAM—MUMBAI
Extending and Modifying LAMMPS
Writing Your Own Source Code
Copyright © 2021 Packt Publishing
All rights reserved. No part of this book may be reproduced, stored in a retrieval system, or transmitted in any form or by any means, without the prior written permission of the publisher, except in the case of brief quotations embedded in critical articles or reviews.
Every effort has been made in the preparation of this book to ensure the accuracy of the information presented. However, the information contained in this book is sold without warranty, either express or implied. Neither the authors, nor Packt Publishing or its dealers and distributors, will be held liable for any damages caused or alleged to have been caused directly or indirectly by this book.
Packt Publishing has endeavored to provide trademark information about all of the companies and products mentioned in this book by the appropriate use of capitals. However, Packt Publishing cannot guarantee the accuracy of this information.
Group Product Manager: Aaron Lazar
Senior Editor: Rohit Singh
Content Development Editor: Rosal Colaco
Technical Editor: Pradeep Sahu
Copy Editor: Safis Editing
Project Coordinator: Deeksha Thakkar
Proofreader: Safis Editing
Indexer: Rekha Nair
Production Designer: Aparna Bhagat
First published: January 2021
Production reference: 1220121
Published by Packt Publishing Ltd.
Livery Place
35 Livery Street
Birmingham
B3 2PB, UK.
ISBN 978-1-80056-226-4
www.packt.com
All praise is to Allah, the most merciful, the most gracious.
– Dr. Shafat Mubin
To all the people who make the world better.
– Jichen Li
Packt.com
Subscribe to our online digital library for full access to over 7,000 books and videos, as well as industry leading tools to help you plan your personal development and advance your career. For more information, please visit our website.
Why subscribe?
Spend less time learning and more time coding with practical eBooks and Videos from over 4,000 industry professionalsImprove your learning with Skill Plans built especially for youGet a free eBook or video every monthFully searchable for easy access to vital informationCopy and paste, print, and bookmark contentDid you know that Packt offers eBook versions of every book published, with PDF and ePub files available? You can upgrade to the eBook version at packt.com and as a print book customer, you are entitled to a discount on the eBook copy. Get in touch with us at [email protected] for more details.
At www.packt.com, you can also read a collection of free technical articles, sign up for a range of free newsletters, and receive exclusive discounts and offers on Packt books and eBooks.
Foreword
Many users of LAMMPS think of it as a tool that provides various features; they likely never look at the source code itself, which is a perfectly fine way to use it. But if you want LAMMPS to do something a little different than what it already does, or you have an idea for something completely new that LAMMPS could do, and you either don't know where to start or want a big-picture view of how to go about it, then this book is for you.
It gives a gentle introduction to how the LAMMPS source code is structured, how to find the code that corresponds to specific features of a LAMMPS input script so you can see how it works, and how to modify or write new code to implement new options or new commands that you want to add to LAMMPS. Most importantly, in my view, it gives some concrete and non-trivial examples of features in different categories (what LAMMPS calls styles) that could be added in the code.
If you follow the step-by-step explanations and instructions provided in this book, you will have a sound understanding of how to do the preceding activities yourself for some new idea you care about in your own research. In fact, about 95% of the LAMMPS source code are styles, which the developers or users added to the code in the same manner and for the same reasons. They wanted to alter an existing model, formulate a new one, or measure a new property, so that LAMMPS could be used to tackle a new problem.
Finally, if you follow the methodologies explained in this book and write some new code for LAMMPS that you think will be useful to others, I hope you will consider submitting it on the LAMMPS GitHub site for possible inclusion in the distribution. That is how many features in current LAMMPS came to be. And it is how its user community continues to help the code grow and become a more useful tool.
Dr. Steven Plimpton
Lead developer of LAMMPS
Contributors
About the authors
Dr. Shafat Mubin (PhD in physics, Penn State) is an assistant professor of physics at Valdosta State University. Since his graduate student days, he has worked with molecular simulations primarily using LAMMPS and has investigated a variety of simulation systems employing a wide array of techniques. He possesses extensive experience in writing custom routines and extending the LAMMPS source code and hosts his own website to instruct and demonstrate how to do so to other users. At present, he is engaged in computational physics research, including molecular simulations, and endeavors to train undergraduate students in computational techniques to help them better prepare for careers in physics.
I want to thank the people who have supported me, especially Farhana.
Jichen Li (graduated from Qingdao University of Science and Technology) is now studying for his master’s degree at the University of Science and Technology of China. He has used LAMMPS to conduct molecular simulations to explore the relationship between polymer microstructure and macro mechanical and rheological properties. He has developed several modeling and post-processing frameworks for LAMMPS and possesses a thorough understanding of its program architecture. He is dedicated to community learning and contributes to the community through his website, where he has written several columns and tutorials for LAMMPS beginners. At present, he is working on trans-scale simulation and the combination of deep learning and simulation.
I would like to thank my parents for their encouragement, Peiyuan Wang, for his support, and my tutor, Jianhui Song, who led me to the field of molecular dynamics.
About the reviewers
Jie Peng is a postdoc from the School of Material Science at Purdue University, US. He has great experience in the area of physics and material engineering. He graduated from the Shanghai Jiao Tong University (China) with a bachelor’s and a master’s degree, and obtained his PhD from the University of Maryland (US) in May 2019. He likes working in interdisciplinary areas where he can see the fundamental science behind materials’ behavior while being able to exploit his knowledge to make use of these materials. He studies the least amount of energy that is transferred in materials – electrons and phonons. The goal is to understand the law of energy transfer at the nanoscale, at which most modern devices operate – cell phones, computers, and so on.
Shih-Hsien Liu is a postdoctoral research associate at the Center for Molecular Biophysics of Oak Ridge National Laboratory. He is a research professional with a chemical engineering background and has experience in physical and computational chemistry for biological applications. His research interests include drug discovery, bioenergy, nanoscience, and molecular dynamics. He obtained his PhD in computational chemistry from the Pennsylvania State University in 2017, an M.S. in materials chemistry from National Taiwan University in 2010, and a B.S.E. in chemical engineering from National Taiwan University in 2008.
Packt is searching for authors like you
If you’re interested in becoming an author for Packt, please visit authors.packtpub.com and apply today. We have worked with thousands of developers and tech professionals, just like you, to help them share their insight with the global tech community. You can make a general application, apply for a specific hot topic that we are recruiting an author for, or submit your own idea.
Table of Contents
Preface
Section 1: Getting Started with LAMMPS
Chapter 1: MD Theory and Simulation Practices
Technical requirements
Introducing MD theory
Understanding the dynamics of point particles
Performing iterative updates using the Velocity Verlet algorithm
Understanding rotational motion
Examining temperature and velocity distribution of particles
Implementing MD simulation practices
Pair-potential cutoff
Periodic boundary conditions
Neighbor lists
Processor communication
Summary
Questions
Chapter 2: LAMMPS Syntax and Source Code Hierarchy
Technical requirements
Introducing the LAMMPS input script structure
Introducing the source code repository
Reviewing the source code hierarchy
Summary
Further reading
Questions
Section 2: Understanding the Source Code Structure
Chapter 3: Source Code Structure and Stages of Execution
Technical requirements
Introducing parent and child classes
fix.cpp and fix.h
pair.cpp and pair.h
compute.cpp and compute.h
Stages of executing the simulation
verlet.cpp
min.cpp
Role of pointers class
Parsing input script commands by input.cpp
Summary
Further reading
Questions
Chapter 4: Accessing Information by Variables, Arrays, and Methods
Technical requirements
Accessing atom properties during simulation runs
Mapping atom indices
Requesting a neighbor list
Accessing physical constants
Reading parameters from the input script
Incorporating new data types
Summary
Questions
Chapter 5: Understanding Pair Styles
Technical requirements
Reviewing the general structure of pair styles
Reviewing the Morse potential
PairMorse::allocate()
PairMorse::settings()
PairMorse::coeff()
PairMorse::init_one()
PairMorse::compute()
PairMorse::single()
Reviewing the table potential
PairTable::settings()
PairTable::coeff()
PairTable::compute()
Reviewing the DPD potential
PairDPD::settings()
PairDPD::compute()
Summary
Questions
Chapter 6: Understanding Computes
Technical requirements
Reviewing the general structure of computes
Reviewing the compute KE class
ComputeKE::ComputeKE()
ComputeKE::init()
ComputeKE::compute_scalar()
Reviewing the compute group/group class
ComputeGroupGroup::ComputeGroupGroup()
ComputeGroupGroup::init() and init_list()
ComputeGroupGroup::pair_contribution()
ComputeGroupGroup::compute_scalar() and compute_vector()
Reviewing the compute RDF class
ComputeRDF::ComputeRDF()
ComputeRDF::init() and init_list()
ComputeRDF::init_norm()
ComputeRDF::compute_array()
Reviewing the compute heat flux class
ComputeHeatFlux::ComputeHeatFlux()
ComputeHeatFlux::init()
ComputeHeatFlux::compute_vector()
Summary
Questions
Chapter 7: Understanding Fixes
Technical requirements
Exploring the general structure of fixes
Reviewing the Fix AddForce class
FixAddForce::FixAddForce()
FixAddForce::init()
FixAddForce::setmask()
FixAddForce::post_force()
Studying the Fix NVE class
FixNVE::initial_integrate()
FixNVE::final_integrate()
Studying the Fix NH class
Studying the Fix Print class
FixPrint::FixPrint()
FixPrint::end_of_step()
Reviewing the Fix Orient/FCC class
FixOrientFCC::FixOrientFCC()
FixOrientFCC::find_best_ref()
FixOrientFCC::post_force()
Analyzing the Fix Wall and Fix Wall/LJ126 classes
Fix wall/lj126
Exploring the Fix Rigid class
FixRigid::initial_integrate()
FixRigid::set_xv()
FixRigid::final_integrate()
Summary
Questions
Chapter 8: Exploring Supporting Classes
Technical requirements
Discovering the Group class
Exploring the Variable class
Studying the Error class
Error::all()
Error::one()
Error::warning() and message()
Reviewing the Memory class
Discussing the Angle and angle/harmonic classes
The Angle class
The angle/harmonic class
Summary
Questions
Section 3: Modifying the Source Code
Chapter 9: Modifying Pair Potentials
Technical requirements
Writing a harmonic potential
Changing class names
Changing variables and equations
Parsing input parameters for Pair Harmonic
Trial run (in.Test_PairHarmonic)
Writing a height-dependent pair potential
Changing the header file, pair_table.h
Parsing input parameters for Pair TableZ
Implementing height-dependence
Trial run (in.Test_PairTableZ)
Writing a friction-based pair style
Modifying the compute() method
Trial run (in.Test_PairCundallStrack)
Summary
Questions
Chapter 10: Modifying Force Applications
Technical requirements
Writing a fix to apply a 2D spring force
Theory (2D spring force)
Modifying the header file – addforceXY.h
Modifying the constructor in fix_addforceXY.cpp
Modifying the init() method in fix_addforceXY.cpp
Modifying the post_force() method in fix_addforceXY.cpp
Trial run (in.Test_FixAddforceXY)
Writing a fix to apply an active force
Theory (active force on bonded atoms)
Modifying the header file – fix_activeforce.h
Modifying the constructor in fix_activeforce.cpp
Modifying the post_force() method in fix_activeforce.cpp
Trial run (in.Test_ActiveForce)
Writing a fix to apply a custom wall force
Theory (expanded Lennard–Jones wall force)
Modifying the header file
Modifying the constructor in fix_wall_region.cpp
Modifying the init() method in fix_wall_region.cpp
Adding the lj126Expanded() method
Modifying the post_force() method
Trial run (in.Test_WallRegion)
Writing a fix to apply a bond-boost potential
Theory (the bond-boost method)
Modifying the header file – fix_addboost.h
Modifying the constructor in fix_addboost.cpp
Modifying the init() method in fix_addboost.cpp
Modifying the post_force() method in fix_addboost.cpp
Trial run (in.Test_FixAddboost)
Summary
Questions
Chapter 11: Modifying Thermostats
Technical requirements
Writing a fix to apply the Andersen thermostat
Theory – Andersen thermostat
Modifying the header file – fix_temp_andersen.h
Modifying the constructor in fix_temp_andersen.cpp
Modifying the end_of_step() method of fix_temp_andersen.cpp
Trial run – in.Test_Andersen
Writing a fix to apply a non-linear temperature increment
Theory – Exponential temperature increment
Modifying the fix_nh.h header file
Modifying the constructor in fix_nh.cpp
Modifying the init() method in fix_nh.cpp
Modifying the compute_temp_target() method in fix_nh.cpp
Modifying the compute_scalar() method in fix_nh.cpp
Trial run – in.Test_Addboost
Writing a fix to print output at evaporation
Theory – Evaporation region
Modifying the header file – fix_printEvaporate.h
Modifying the constructor and init() methods in fix_printEvaporate.cpp
Modifying the end_of_step() method in fix_printEvaporate.cpp
Trial run – in.Test_PrintEvaporate
Summary
Questions
Appendix A: Building LAMMPS with CMake
Technical requirements
Prerequisites for working with LAMMPS
Downloading the source code
Installing the dependencies
Downloading MPICH
Building LAMMPS
Including packages in the build
Including modified codes
Building with CMake
Compiling LAMMPS in Windows
Downloading Git for Windows
Downloading the source code
Downloading CMake for Windows
Downloading Microsoft MPI
Downloading Code::Blocks
Developing LAMMPS
Building with Make
Further reading
Appendix B: Debugging Programs
Technical requirements
What is debugging?
Prerequisites
Method 1 – Debugging with GDB
Method 2 – Debugging with VSCode
Insight into sbmask()
Further reading
Appendix C: Getting Familiar with MPI
Technical requirements
What is MPI?
MPI initialization
MPI finalization
MPI current process ID
MPI (number of processes)
MPI send message
MPI receive message
MPI message
MPI in LAMMPS
Example of evaluating
Data exchange between owned atoms and ghost atoms
Appendix D: Compatibility with Version 29Oct20
Technical requirements
Translating Fix Widom into the 3Mar20 version
Incompatibility 1 – Parsing input parameters
Incompatibility 2 – Atom class
Incompatibility 3 – Function signatures in group.cpp
Assessments
Other Books You May Enjoy
Preface
LAMMPS is one of the most widely used tools for running simulations for research in molecular dynamics. While the tool itself is fairly easy to use, more often than not, you'll need to customize it to meet your specific simulation requirements.
Extending and Modifying LAMMPS bridges this learning gap and helps you achieve this by writing custom code to add new features to LAMMPS source code. Written by ardent supporters of LAMMPS, this practical guide will enable you to extend the capabilities of LAMMPS with the help of step-by-step explanations of essential concepts, practical examples, and self-assessment questions.
Filled with real-world examples, this book teaches you how to find your way around the LAMMPS source code and shows you, in a pragmatic and easily digestible way, how to modify the code to achieve custom simulations for your work.
This LAMMPS book provides a hands-on approach to implementing associated methodologies that will get you up and running and productive in no time. You'll begin with a short introduction to the internal mechanisms of LAMMPS and gradually transition to an overview of the source code along with a tutorial on modifying it.
You'll also learn how to identify how LAMMPS input script commands are executed within the source code and understand the architecture of the source code. You'll learn how to modify the source code to implement custom features in LAMMPS, relate source code elements to simulated quantities, and explore pair styles and computes within LAMMPS. You'll also learn how to create various custom fixes to control your simulation using pre-defined criteria.
As you advance, you'll learn about the structure, syntax, and organization of LAMMPS source code, and will be able to write your own source code extensions to LAMMPS that implement features beyond the ones available in standard downloadable versions. You'll understand how to work with custom solutions for applications on force, thermostats, and pair potentials.
By the end of this source code book, you'll have learned how to how to add your own extensions and modifications to the LAMMPS source code that can implement features that suit your simulation requirements.
Who this book is for
This book is for students, faculty members, and researchers who are currently using LAMMPS or considering switching to LAMMPS, have a basic knowledge of how to use LAMMPS, and are looking to extend LAMMPS source code for research purposes. This book is not a tutorial on using LAMMPS or writing LAMMPS scripts, and it is assumed that the reader is comfortable with basic LAMMPS syntax. The book is geared toward users with little to no experience in source code editing. Familiarity with C++ programming is helpful but not necessary.
What this book covers
Chapter 1, MD Theory and Simulation Practices, briefly introduces fundamental concepts of molecular dynamics, common techniques, and features that will be addressed later in the book to provide a link between theory and simulation.
Chapter 2, LAMMPS Syntax and Source Code Repository, introduces LAMMPS input script and the execution of MD simulation features through LAMMPS commands and an introduction to the repository of standard LAMMPS code.
Chapter 3, Source Code Structure and Stages of Execution, introduces the LAMMPS source code format and the sequence of execution followed in a LAMMPS script.
Chapter 4, Accessing Information by Variables, Arrays, and Methods, describes the elements that store information and perform calculations during a simulation run.
Chapter 5, Understanding Pair Styles, describes and dissects selected pair styles.
Chapter 6, Understanding Computes, describes and dissects selected computes.
Chapter 7, Understanding Fixes, describes and dissects selected fixes.
Chapter 8, Exploring Supporting Classes, describes and dissects other additional LAMMPS features.
Chapter 9, Modifying Pair Potentials, discusses the implementation of custom pair styles.
Chapter 10, Modifying Force Applications, discusses the implementation of custom fix styles to apply forces based on pre-defined criteria.
Chapter 11, Modifying Thermostats, discusses the implementation of custom fix styles associated with temperature control.
Appendix A, Building LAMMPS with CMake, introduces LAMMPS compilation using CMake in Linux and Windows.
Appendix B, Debugging with GDB and Visual Studio Code, introduces debugging tools with applications on understanding bit operations.
Appendix C, Getting Familiar with MPI, introduces MPI parallel processing.
Appendix D, Compatibility with Version 29Oct20, discusses source code differences between stable versions 3Mar20 and 29Oct20.
Assessments, contains the answers to all the questions from this book.
To get the most out of this book
If you are using the digital version of this book, we advise you to type the code yourself or access the code via the GitHub repository (link available in the next section). Doing so will help you avoid any potential errors related to the copying and pasting of code.
Also, we advise that you download the LAMMPS stable version 3Mar20 and have the corresponding files open when viewing the code discussions provided in the book.
Download the example code files
You can download the example code files for this book from GitHub at https://github.com/PacktPublishing/Extending-and-Modifying-LAMMPS-Writing-Your-Own-Source-Code. In case there's an update to the code, it will be updated on the existing GitHub repository.
We also have other code bundles from our rich catalog of books and videos available at https://github.com/PacktPublishing/. Check them out!
Download the color images
We also provide a PDF file that has color images of the screenshots/diagrams used in this book. You can download it here: https://static.packt-cdn.com/downloads/9781800562264_ColorImages.pdf.
Conventions used
There are a number of text conventions used throughout this book.
Code in text: Indicates code words in text, database table names, folder names, filenames, file extensions, pathnames, dummy URLs, user input, and Twitter handles. Here is an example: "In the single() method, the scalingZ variable is calculated and used to scale the potential before it is returned."
A block of code is set as follows:
units metal
dimension 3
boundary p p p
atom_style atomic
atom_modify map array
When we wish to draw your attention to a particular part of a code block, the relevant lines or items are set in bold:
pair_style tableZ linear 10000
pair_coeff 1 1 Table_Potential.table V_Table 12.0 17.0 1.0 -0.04 10.0
compute TABLE_PE G1 group/group G2
variable TABLE_PE equal c_TABLE_PE
Bold: Indicates a new term, an important word, or words that you see onscreen. For example, words in menus or dialog boxes appear in the text like this. Here is an example: "As you can see, the atoms will cross R_EVAP sequentially and the Fix PrintEvaporate command will be able to register these crossings by printing to file."
Tips or important notes
Appear like this.
Get in touch
Feedback from our readers is always welcome.
General feedback: If you have questions about any aspect of this book, mention the book title in the subject of your message and email us at [email protected].
Errata: Although we have taken every care to ensure the accuracy of our content, mistakes do happen. If you have found a mistake in this book, we would be grateful if you would report this to us. Please visit www.packtpub.com/support/errata, selecting your book, clicking on the Errata Submission Form link, and entering the details.
Piracy: If you come across any illegal copies of our works in any form on the Internet, we would be grateful if you would provide us with the location address or website name. Please contact us at [email protected] with a link to the material.
If you are interested in becoming an author: If there is a topic that you have expertise in and you are interested in either writing or contributing to a book, please visit authors.packtpub.com.
Reviews
Please leave a review. Once you have read and used this book, why not leave a review on the site that you purchased it from? Potential readers can then see and use your unbiased opinion to make purchase decisions, we at Packt can understand what you think about our products, and our authors can see your feedback on their book. Thank you!
For more information about Packt, please visit packt.com.
Section 1: Getting Started with LAMMPS
A LAMMPS input script performs molecular dynamics simulations according to the input commands provided. But what happens behind the scenes as each command is parsed and executed? How are these commands processed and how are relevant calculations invoked within the source code?
This section introduces you to molecular dynamics theory and standard simulation practices that are implemented in LAMMPS, along with the corresponding syntax. The theory and practices will be revisited in later sections when the LAMMPS source code is analyzed.
This section comprises the following chapters:
Chapter 1, MD Theory and Simulation PracticesChapter 2, LAMMPS Syntax and Source Code RepositoryChapter 1: MD Theory and Simulation Practices
This chapter introduces the theory behind molecular dynamics (MD) and some common simulation practices. Starting from Newton's laws, we outline the physics behind the dynamics of point particles and rigid bodies, discuss iterative updating and the relevance of temperature, and end by listing computational practices.
In this chapter, we will cover the follow topics:
Introducing MD theoryUnderstanding the dynamics of point particlesPerforming iterative updates using the Velocity Verlet algorithmUnderstanding rotational motionExamining temperature and velocity distribution of particlesImplementing MD simulation practices including cutoff, periodic boundaries, and neighbor listsBy the end of this chapter, you will have grasped an understanding of the theoretical fundamentals implemented in MD software.
Technical requirements
You can find the full source code used in this book here: https://github.com/PacktPublishing/Extending-and-Modifying-LAMMPS-Writing-Your-Own-Source-Code
Introducing MD theory
MD is based on simulating individual particle trajectories over a desired time period to analyze the time evolution of an entire system of particles in a solid, liquid, or gaseous state. Each particle (usually an atom) is allowed to traverse in space as determined by Newton's laws of classical dynamics, where the atomic positions, velocities, and accelerations at one point in time are used to calculate the corresponding kinematics quantities at a different point in time. This process is repeated over a sufficiently long-time interval for every atom in a system, and the final configuration of the atoms indicates the time evolution of the system over the said time interval.
Typical MD simulations are limited to the study of atomistic systems consisting of atoms in the range of to , occupying a simulation box with a length in the order of nanometers, over a regular timescale of nanoseconds. MD simulations of such microscopic systems are relevant when the systems are able to represent the time evolution of corresponding macroscopic systems.
The theory behind MD is described briefly in the following sections. For a more detailed understanding, you are advised to refer to the dedicated literature on MD theory (Computer Simulation of Liquids by Michael P. Allen and Dominic J. Tildesley, and Understanding Molecular Simulation by Daan Frenkel and Berend Smit).
Understanding the dynamics of point particles
The trajectory of a point particle over time t is calculated from its mass m and net force using Newton's equation, as illustrated here:
is determined from the sum of all forces acting on the particle. In a system of interacting particles, the force acting between one pair of particles can be determined from the gradient of the potential energy function . Here, is the displacement vector that points to the particle for which we are calculating (using the following equation) and originates from the other particle in the pair (shown in Figure 1.1):
This gives us the three components of . The (x) force component is given as follows:
The (y) force component is calculated by the following formula:
The (z) force component is given by the following formula:
Here, is the distance between the pair of particles. By Newton's third law, it follows that the force components acting on the other particle in the pair have the same magnitudes, with the opposite sign.
The following diagram illustrates this concept, using two particles interacting via a 12-6 Lennard-Jones potential of the form , where represents the well depth and represents the position at which the potential is zero:
Figure 1.1 – Two particles in 1D interact via a Lennard-Jones potential
In this diagram, if we want to calculate the force from particle 2 (located at ) acting on particle 1 (located at ), then we use . Subsequently, the reaction force acting on particle 2 is given by .
Since potential energy functions are commonly expressed as functions of r, the expressions we saw earlier make it more convenient to calculate the force components. The sum of forces acting on a single particle by its interaction with all particles in the system gives the net force, as illustrated here:
Altogether, we obtain three equations to solve for , and . The x(t) equation is given as follows:
The y(t) equation is given by the following formula:
The z(t) equation is given by the following formula:
These values are used to generate the complete trajectory of the particle over a desired time interval. This process is repeated for all particles in the system to yield the complete system time evolution in the same time interval.
Performing iterative updates using the Velocity Verlet algorithm
As particles move in a system over time, the separation between each pair of particles changes accordingly, resulting in a change of the derivatives described in the previous section—that is, the force becomes time-dependent. Therefore, the trajectory is updated over a small increment of time called the timestep, and by iteratively updating the trajectory over a large number of timesteps, a complete trajectory can be obtained.
The purpose of keeping the timestep small is to ensure that the particles do not undergo drastic changes in position and therefore conserve energy in the system. For atomic masses, especially in the presence of forces from molecular bonds, angles, and dihedrals, a timestep of approximately 1 femtosecond is typically employed.
An advantage of using a small timestep is that the net force acting on a particle can be approximated to remain constant in the duration of the timestep. Therefore, the equations of motion that iteratively update the trajectory at timestep increments become the following:
Here, represents the velocity vector of the particle, and its components are iteratively updated by the following:
Here, the terms within the [ ] represent the average acceleration in each dimension, calculated using the accelerations at time t and the following iteration, . This is known as the Velocity Verlet algorithm, which considerably reduces errors over a long simulation period as compared to algorithms that use the acceleration at a single point in time (for example, the Euler algorithm). Furthermore, the Velocity Verlet algorithm is able to conserve energy and momentum within rounding errors with a sufficiently small timestep, unlike the Euler algorithm, which can lead to an indefinite increase in energy over time.
In effect, the position and velocity of each particle are tabulated at each iteration, as illustrated in the following table:
Table 1.1 – Table showing sequence of iterative update of position and velocity vectors of each particle
This table shows the sequence of iterative update of point particles that undergo a linear motion without any rotational component. In the case of a non-point object such as a rigid body, both linear and rotational motion must be incorporated for a proper treatment of its dynamics. Similar to linear motion, rotational motion of rigid bodies can be iteratively updated using a similar algorithm, as discussed next.
Understanding rotational motion
Rigid bodies are often used to represent molecules with inactive vibrational degrees of freedom. In a rigid body, the locations of all constituent atoms are frozen with respect to the rigid-body coordinates. The intramolecular forces from bonds, angles, dihedrals, and impropers are assumed not to create any distortion of the rigid body. Therefore, for the purposes of MD, all intramolecular forces in a rigid body can be disregarded. Any force exerted on any constituent atom in a rigid body acts on the entire rigid body. In addition to changing its linear velocity, this force can create a torque on the rigid body and change its angular momentum and angular velocity. The total force on a rigid body composed of N particles is calculated using the following formula:
Using the displacement vector of particle i measured from the center of mass of the rigid body as the origin, the torque of the rigid body is calculated by the following formula:
The angular momentum of the rigid body is obtained from the time-integration of using the Velocity Verlet algorithm, as follows:
The angular velocity of the rigid body is obtained by a two-step procedure from using an intermediate step at a half-timestep and a second step at a full-timestep, which will be described in more detail when analyzing the dynamics of rigid bodies in Chapter 7, Understanding Fixes. The procedure is illustrated here:
Here, is the moment of inertia tensor of the rigid body.
For a system of point particles or rigid bodies, the average kinetic energy and the average speed or angular speed are determined by the system temperature. Therefore, controlling the temperature adequately is often essential when simulating a molecular system. The next section discusses the features of a molecular simulation determined by the temperature.
Examining temperature and velocity distribution of particles
A system at thermal equilibrium at a constant temperature T is characterized by its Maxwell-Boltzmann velocity distribution. According to this distribution, the probability distribution of velocities in a single direction i (which can be x, y, z) of a system of particles of mass m each is given by the Gaussian function, illustrated here in Figure 1.2:
Figure 1.2 – The Maxwell-Boltzmann velocity distributions (left) and speed distributions (right)
The preceding graph is plotted for the same system at three different temperatures .
The corresponding functional form that depends on the mass and temperature is shown here:
Here, is the Boltzmann constant. This distribution has a mean of and a standard deviation of . The shape of the Gaussian curve is determined by the ratio of . The velocity distribution of the velocity vector is given by the following formula:
In spherical coordinates, this distribution can be written in terms of the speed , as follows:
This is the Maxwell-Boltzmann speed distribution, also known as a Rayleigh distribution. The shape of the speed distribution changes with temperature, as shown in Figure 1.2, and the peak speed increases with temperature. The velocity distributions are wider at higher temperatures, and the speed distributions show larger peak speeds at higher temperatures. An algorithm that controls temperature in a molecular simulation must account for the preceding features regarding particle velocities.
So far, concepts that dictate the operation of an MD simulation have been discussed. These concepts are implemented in a computational environment through codes that will be discussed in later chapters. In the next section, we will discuss related computational concepts that are commonly encountered in MD simulations and that are prevalently used to enhance simulation performance.
Implementing MD simulation practices
In order to implement MD simulations that are computationally efficient and model realistic atomic systems, a number of standard simulation practices are employed. A user is likely to find these practices employed in typical Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) scripts and it is therefore helpful to be familiar with these concepts before delving into the LAMMPS source code. In this section, we will briefly discuss some of these practices (details are available in the MD textbooks referred to in the Introducing MD theory section).
Pair-potential cutoff
If it is desired that a potential reaches exactly zero at the cutoff, an offset () can be employed by calculating the potential at the cutoff , that is . This offset can then be subtracted from the original potential to guarantee a zero value at the cutoff, that is . Altogether, the potential with the offset changes the value of the system potential, but does not alter the forces because is a constant term that produces zero force contribution upon differentiating.
Most pair-potential functions are defined to asymptotically approach zero potential— for example, the Lennard-Jones function approaches zero as the inverse-sixth power of separation. For particles located far from each other, it implies that there exists some small but non-zero potential between them, which may be negligible in a simulation context but still adds numerical computation overhead.
Therefore, a common practice is to employ a cutoff distance for pair potentials, beyond which the interaction is assumed to be zero. The cutoff is chosen at a separation where the potential value is sufficiently small so as not to create significant errors by discounting interaction with neighbors located beyond the cutoff. During a simulation run, only the neighbors located within the cutoff radius of a particle are considered for force or potential calculations, thereby reducing computation time.
Periodic boundary conditions
Simulated systems often consist of a small cell that is representative of a larger system—for example, nanocrystal representation of a metal lattice. In such systems, it can often be assumed that the entire system is made up of many replications of the simulation box. In 1912, Born and von Karman came up with the implementation of periodic boundary conditions that can simulate a continuous replicated system starting from a simulation box.
Each wall (that is, boundary) of the simulation box is assumed to be adjacent to an identical simulation box containing identical atoms at identical positions. In effect, each wall of the box leads to a replica of the simulation box, containing images of the same atoms as in the simulation box. At every iteration, as the atom positions in the simulation box are updated, the image atoms in the replicas are accordingly updated. The following diagram shows this simulation box:
Figure 1.3 – A 2D simulation box without (top) and with (bottom) periodic boundaries
In this diagram, we see the following:
Top: A snapshot of a 2D simulation box consisting of four particles (indexed from 1 to 4) without periodic boundaries.Bottom: The same simulation box with periodic boundaries in all directions produces 8 replicas (26 replicas in 3D) containing image atoms (dashed circles). At every iteration, all replicas are identically updated by replicating the particles in the central simulation box.This way, enclosing wall boundaries are eliminated, and any atom that exits the simulation box by passing through a wall will automatically re-enter through the opposite wall, thus conserving the number of atoms and modeling continuous motion of atoms. Furthermore, the atoms in the simulation box are permitted to interact with the image atoms located in the replicas, ensuring that atoms located near to the walls are not considered as fringe atoms with few neighbors only because of their remoteness from the box center.
In Figure 1.3, in the absence of periodic boundaries, particle 3 at the top-right corner is observed to be located far from particles 1 and 4, and would undergo reduced or zero interactions with these particles. By contrast, when periodic boundaries are implemented, particle 3 is located considerably closer to images of particles 1 and 4.
Subsequently, particle 3 is able to interact with other particles located in all directions around itself as if it were contained in a continuous, unbounded simulation box. At the same time, consistency of interaction requires that particle 3 interacts with either an image or with the original particle in the simulation box, but not both together, and it is accomplished by setting the side lengths of the simulation box to be at least double the pair-potential cutoff distance between atoms. This requirement is known as the minimum image convention, which guarantees that pair-potential interactions are not double-counted.
Long-range electrostatic potentials decay considerably slower than many other interactions (as inverse of the separation) and would therefore require disproportionately long cutoffs and correspondingly large simulation boxes to model with a traditional pair potential. To circumvent this problem, electrostatic interactions with image atoms are summed up by an Ewald summation (or a related particle-mesh method) that divides the calculation to a real-space component and a reciprocal-space component. This way, a cutoff is not required, but periodic boundaries are necessary to ensure sum convergence.
When an atom exits a simulation box and its image enters from the opposite side of the box, the atom coordinates can extend beyond the simulation box coordinates. This is accounted for using the concept of wrapped and unwrapped coordinates, where unwrapped coordinates represent the unadjusted atom coordinates and wrapped coordinates represent the coordinates adjusted by resetting to the coordinates of the re-entry wall.
In LAMMPS, trajectory output files may include an image flag to keep track of wrapped and unwrapped coordinates. When an atom exits the simulation box along the positive direction in any dimension, the image flag corresponding to this axis would increment by one. Accordingly, the image flag decrements by one if the atom exits along the negative direction. Thus, the image flag of an atom multiplied by the corresponding simulation box side length can be used to convert from wrapped to unwrapped coordinates.
Neighbor lists
An atom only interacts with its neighbors located within a cutoff radius, and these neighbors are identified by calculating their distances from the atom in consideration. For a system of N atoms, this could lead to calculating the distance between each of the pairs of atoms. To reduce the computation overhead, for every atom a subset of neighbors is selected into a neighbor list (suggested by Verlet in 1967), and only these short-listed neighbors are used to calculate the interactions with the atom.
At the beginning of a simulation, a neighbor list is built for every interacting atom in the system by tagging all its neighboring atoms that lie within its cutoff or within a short buffer width, known as the skin width (shown in Figure 1.4). If the atoms are not traveling at extreme speeds, only the atoms located within the cutoff or the skin (that is, the neighbor-list members) at a certain iteration can be expected to lie inside the cutoff radius in the next iteration, and it can be expected that no atom previously located outside the skin will be able to cross inside the cutoff radius in the space of one iteration.
This way, the neighbor list excludes atoms that are located sufficiently far by reusing information from the preceding iteration, and this therefore reduces computation time in calculating neighbor distances. The process is illustrated in the following diagram:
Figure 1.4 – Illustration of the radial cutoff distance and the skin width of a central atom (black dot)
At any timestep, the atoms that are located inside are included in the neighbor list of the central atom, whereas only the atoms located inside interact with the central atom. At the next iteration, only the atoms tagged in the neighbor list are considered when identifying atoms that can interact with the central atom, and the neighbor list may be rebuilt depending on the displacements of the atoms.
Generally, the neighbor lists do not have to be rebuilt at every iteration, and the interval between rebuilding can be defined by the user. Also, a common practice is to rebuild only when an atom has traveled more than half the skin width since the previous neighbor-list build. It also implies that a larger skin width requires less frequent neighbor-list builds, at the cost of including a greater number of neighbors per list at every iteration.
The skin width can be specified independently of the cutoff distance and is generally chosen according to timestep and expected particle velocities. At the lower limit, the skin width should be double the maximum displacement expected of a particle in an iteration. In LAMMPS, if an atom is able to cross the skin width in one iteration, the associated neighbor lists experience a dangerous build whereby the loss of the atom can lead to errors in force calculations and violation of energy conservation. If a dangerous build is detected, the neighbor list needs to be rebuilt to rectify this.
When a neighbor list is used to calculate the interaction force between a pair of atoms, the two atoms in the pair can be assigned equal and opposite forces by Newton's third law. This equivalency can be enabled using a full-neighbor list, whereas it can be disabled by using a half-neighbor list. While a full-neighbor list reduces computation cost, half-neighbor lists may be preferred in simulations where equal and opposite forces are not applicable.
Processor communication
Modern parallel computers have two main forms: shared-memory machines, where multiple processors (often called cores) access the same memory; and distributed-memory machines, where a processor (often called a node) cannot access the memory designated to another processor. Modern high-performance computing facilities are all based on a hybrid architecture. A node consists of multiple cores, and many nodes combine to form a supercomputer. This architecture also leads to an issue of distributing the workload and allocating tasks between nodes and cores.
In the context of MD, one strategy that addresses the task allocation problem is spatial decomposition, which divides the simulation box into equal sections and assigns atoms located in each section to a different core, as shown here:
Figure 1.5 – (Left) A single processing core calculating the trajectory of every atom in the simulation box. (Right) The simulation box is divided into domains, and the atoms in each domain are assigned to a different core
The domains in the preceding diagram are spatially decomposed by volume. Copies of atoms, known as ghost atoms, are exchanged between processors to account for interactions between atoms located on different cores. Each core calculates the interactions of atoms assigned to it in parallel and thereby increases the overall simulation speed.
To account for interactions between atoms located on different cores, each core builds its domain with a shell that accommodates particles located at the edges of a domain, and exchanges copies of particles (known as ghost atoms) with other cores, as illustrated in the following diagram:
Figure 1.6 – A core domain showing its shell (dashed box) and ghost atoms shared with a neighboring domain
There are two stages of communication between cores involved, detailed here:
In the first stage, the shell exchanges the ghost-atom coordinates with neighboring cores to register the particles that can interact with the particles contained in the domain. The ghost atoms in these shells are essentially images of atoms in other domains.In the second stage, the updated atom positions at the following iteration are used to determine whether any atom has moved to a different domain, and if so, the atom coordinates and kinematics information are communicated to the appropriate cores. After this step, the particles in the shell are deleted from the memory.This parallel mechanism demonstrates that communication time between cores multiplies as the number of cores increases. Due to memory and bandwidth limits, parallelization cannot achieve the ideal optimization of , where N is the number of atoms and P is the number of processors, and leads to a sub-linear speedup. Therefore, for a given simulation system an increasing number of cores eventually reduces efficiency, and there exists an optimum number of cores that delivers the best performance.
Summary
The basics of MD and common MD simulation practices have been laid out in this chapter to help elucidate the physics and mathematics implemented in MD codes. When the LAMMPS source code is discussed in later chapters, the concepts discussed here will be referenced and explicated in terms of the source code.
In the next chapter, you will be introduced to LAMMPS input script and the execution of MD simulation features through LAMMPS commands, and the repository of standard LAMMPS codes.
