Tutorial | LAMMPS 1.0

LAMMPS 入门

LAMMPS为理论与计算化学常用软件代表，类似的还有VASP， Gaussian，CASTEP, Gromacs等等，这里仅以LAMMPS作为入门举例。由于科研内容的复杂度高、计算量大，这类软件有几个特点：1）一般都需要在linux系统下的服务器上多核运行 2）有输入文件和输出文件，中间为黑箱，基本结构和语法相似。其中产生输入文件和阅读输出文件都需要其他软件或程序，需要一起学习。3）主要学习手段为软件官方用户手册，里面对所有命令有使用方法、相关物理原理、引用文献。这是和其他training中自编软件的主要区别：自己写程序只需要基本物理原理，其他自己都是清楚的、自定义的；而程序需要你熟悉他人程序语言框架，但后期更省时间

step 0 : linux环境熟悉+分子动力学基础

  • 安装winscp后， 安装putty，地址114.214.200.31 ，登录学习账号xue, 密码找我要（后期会给你开单独账号），此软件使用与window文件管理类似。自建属于你的文件夹，之后操作都在自己文件夹下，不要用别人的。点击putty图标，再次登录，此为shell会话。在shell会话下，熟悉新建文件夹、复制、删除文件、浏览文件和修改文件内容相关命令；

    ref: https://blog.csdn.net/…./129345242 第（二）部分。需要记住3-7，其他只要跟着尝试一下即可，即用即查

    也可以随便B站CSDN搜Linux入门

  • 分子动力学基本原理组LAMMPS。分子动力学基本原理和方法可以参考training-活性部分。或学习以下视频：

    ref： B站的lammps入门介绍

    学习完成后回答：

    • 什么是周期性边界？
    • 什么是NVT, NPT, NVE系综？
    • 分子动力学的牛顿方程怎么写，每一项在LAMMPS中存储在哪里？
step 1 : 体验第一个lammps程序

目标： 看到一个10条链的动图。以我个人理解，

  • /test/01下有个简单的程序，打开shell 会话，进入该目录，提交命令为 qsub sub.pbs 提交后用qstat查看状态，有个R就是在运行了，c就是运行结束（也可能是程序死了）。
    • 为了防止误删,建议养成习惯，每次跑新任务先复制新建文件夹，做好编号。每次提交任务，都要新建一个文件夹并编号，以免报错忘了自己改过什么；及时清理dump出来的巨大文件，可以保留输入文件，即使程序报错了，以便后面参考， 误删了也可以重跑。类似于做实验写实验记录， 所有成功失败都要可以参考。
  • 如果跑成功了会生成 .xyz文件，导入OVITO软件即可查看动图,可保存为gif；运行过程随时可打开. log 或log.lammps文件查看进度，如果程序死了也是在这两个文件结尾看报错。如果都没有报错再看err
  • 保持良好习惯，做实验记录：计算机实验也需要记录，以文件夹标号，简单记录跑了什么体系，有哪些参数，看到了什么。可用word，text或markdown(极力推荐，科研利器)
  • 学习完成后回答：
    • 成功提交一个LAMMPS任务需要哪些文件？分别包含哪些部分？
    • in文件中的每一行是干什么的？（step 0课程中有部分答案）
    • 注意 pair_style lj/cut的参数， 1.0 1.0 1.12代表什么？

(ref: https://docs.lammps.org/atom_style.html)，其他命令也可在手册中搜索达到

  • 新建文件夹后，1) 修改为周期性边界2) 将链内珠子的相互作用从纯排斥改为吸引（提示只需要修改一个参数），再次运行，查看结果、记录实验。
step 2: 跑两条真实无序蛋白的粗粒化链

目标： 参考test用c++或者python等，写个小程序，能够通过读入一个序列输出datafile 中原子坐标， 和连键的信息，能够运行

编程是重要的科研工具，无论你做哪个方向都是绕不开的，此处只需要针对问题进行熟悉代码（遇到什么不会查什么），不必从头学某个语言。注意，这个练习只是为了让你适应一点编程思维，不是为了让你写出多好用多高效的程序。一共就几十行即可完成，不要查网络上实现这个功能的完整代码，自己拿空文本写然后遇到问题问我或者查那个问题就行

  • 从数据库下载某个感兴趣的无序蛋白的序列保存至txt文件（如 uniprot.org是常用蛋白数据库）。

  • 写一个生成datafile的小程序，在程序中实现将序列读入数组、并判断每个氨基酸的带电。datafile基本产生思路：

    • （0）开二维数组，保证能存储所有粒子的坐标

    • （1）每条链的第一个原子用随机数确定坐标

    • （2）后面的每个原子生成一个随机数确定它和上一个原子的矢量（如果两个原子相切，三维距离应该是二者半径之和）

    • （3）需要判断是否已经生成的原子中有和新原子距离过近的（可以相互接触，不要距离是0.0几这种就行，容易导致程序发散），要重新生成

    • （4）最后对于生成原子超出盒子边界的，需要根据周期性边界的原理把它约化到盒子内

    • （5）仅在输出时用到粒子类型、电荷的数组，根据氨基酸的电荷区分粒子类型 1-不带电 2-正电 3-负电

    • 一些减少你工作量的代码：

      1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      20
      21
      22
      23
      24
      25
      26
      #include <fstream>
      #include <iostream>

      fp2 = fopen(datafile_name, "w");

      fprintf(fp2,"LAMMPS data file for HSF. generated by c++\n\n");
      fprintf(fp2, "       %d   atoms\n", totalp*bead+ totalcation);
      fprintf(fp2, "       %d   bonds\n\n", totalp * (bead - 1));
      fprintf(fp2, "          21   atom types\n");
      fprintf(fp2, "          1   bond types\n");
      fprintf(fp2, "  %f       %f        xlo xhi\n", 0.0f, box_len);
      fprintf(fp2, "  %f       %f        ylo yhi\n", 0.0f, box_len);
      fprintf(fp2, "  %f       %f        zlo zhi\n\n", 0.0f, box_len);
      fclose(fp2);

      fscanf(fp1, "%c", &atomtype[i]);

      double rn()
      {
      	int x, y;
      	double z;
      	x = rand();
      	y = x % 1000;
      	z = y / 1000.0;
      	return(z);  
      }

  • data文件还需要（1）bond模块 （2）mass模块， 参考step1 的举例完成

  • 注意到我们需要让datafile包含电荷信息，这在step1里的datafile中是没有的。因此需要把in文件的atom_style 改成real, 并修改相应的datafile格式

    ref :https://docs.lammps.org/atom_style.html (有点复杂，不如百度)

  • 此处我们将采用residue-based 粗粒化模型，每个氨基酸用一个球模拟。每种氨基酸球有独立的电荷（参与电荷相互作用）、疏水参数lambda(参与LJ短程势)。

    ref： https://doi.org/10.1016/bs.mie.2020.07.009

    因此需要将in文件pair style 改为lj/cut/coul/debye. 思考怎么改参数。可百度lammps in文件每一列的用途，和我给的test唯一区别是原来只有一种粒子“1”，现在有三种“1，2，3”。了解命令dielectric。

    ref:pair_style lj/cut/coul/debye

  • 一定会报错，这很正常，尝试自己分析原因，然后问我

  • 学习完成后回答：

    • 蛋白在生理盐水中。溶剂/离子体现在哪里？
    • 了解lj/cut/coul/long 和lj/cut/coul/debye. 注意到我之前用的是long, 但我计划让你用debye, 思考下为什么我们选用debye（包括科学性和可行性）
    • 电荷作用/lj的截断应该设在半径的多少倍多少？为什么是截断？
    • 体系为什么需要warm up? 有哪些warm方法？（如果你遇到了lost atom、bond/atom 报错，那就很可能是因为没有warm up，思考为什么lost atoms）

2条成功后，可以尝试跑个100条的。成功后，恭喜你，此后为输入输出省下大量时间。

step 3 : 蛋白的相分离与统计分析

目标： 熟悉lammps其他in文命令（重点熟悉fix)，熟悉一些统计内容和方法（是的，又要写统计的程序了）。

  • 程序将变得复杂。相对于step2 里的in文件，你需要学习的内容（可以查手册了解详情）：
    • variable 命令是什么，怎么替换变量的（$符号是linux系统中的变量替换符号）；loop命令实现循环的方法（类似于for循环）
    • lammps的运行顺序是按in文件逐行读取，因此可以出现多个run, 执行某个run时后面的命令不会被读到。观察理解为什么需要多个run
    • 再次熟悉dump命令, 里面有控制跑多少步次输出一次的， 设置大一些，要保证你run完整个xyz文件不要超过500帧， 一开始甚至100帧就足够了。包含所有原子坐标的文件极大，不能浪费空间。
    • 真正要跑相分离，需要更多的步数，可以考虑再加一两个0

如果你理解了以上，则能看懂下面的命令，并添加：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
。。。。。。
run 10000000#(step 2 的最后)

#循环升温
label loopaaa# 循环入口

variable i loop 10 
variable t equal $t+0.1
 
fix		10 all nvt temp $t $t 100.0
dump	$i all atom ${frame1} temp_$i.dump
run 5000
undump  $i 
unfix 10

next i
jump hsf.lmp loopaaa

尝试让它跑完，观察输出。

  • 理解direct coexistence 方法和slab构型的意义。

    ref: 主要参考文献还是上面Mittal的那篇，第二章疏水不用看，主要是第三章的slab method。大佬写的教学版， 非常通俗易懂，里面包括lammps里怎么构建slab模型， 和后面跑完之后怎么数据处理。

    接下来，lammps中实现slab。 基本思路是： 原来的datafile即可, minimize→平衡→压缩xy→拉长z→平衡→进入升温循环。升温循环已写，改变盒子的代码如下：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    #1. 让蛋白聚集
    fix           1 all nvt temp  $t  $t 100.0 
    run             1000
    unfix 1

    #2. 将盒子沿着zy方向压缩
    fix		3 all nvt temp $t $t 100.0
    fix     4 all deform 10 x final 0.0  150.0  y final 0.0  150.0 
    run 1000
    unfix 3
    unfix 4

    #3.拉长时间
    【手动打码】和2的压缩写法一模一样，沿着z拉伸。
    run 6000
    unfix 7
    unfix 8

    请添加至in文件合适的部分

    如果运行成功，你已经实现了模拟的全部步骤。行百里路___90%了，加油！

  • 最后，统计最大cluster, 得到相图。理论上，只要有dump文件里的粒子信息（以及速度信息，这里只是没有输出），我们就知道了体系的全部微观信息，从这里可以统计任何宏观力学量。这显然并不容易，你要写很多代码，好在lammps内置了很多统计程序，可以省很多麻烦。你需要学习：

  • compute chunk/atom, 理解chunk,

  • fix命令下有很多子命令，是改变体系性质，或者改变输出内容，了解这里新出现的fix ave/time。

1
2
3
4
5
6
7
8
9
10
11
12
13
#添加统计量
#-----------------------density distribution-----------------------
compute 	den1 all chunk/atom bin/3d x lower 3  y lower 210   z lower 210 
compute 	myChunk1 all property/chunk  den1 count
#max cluster
compute 	cluster1 par cluster/atom 7.0
compute 	cc1 par chunk/atom c_cluster1 compress yes
compute 	size par property/chunk cc1 count


#输出统计量
fix clust1 all ave/time 1 1 5000 c_size c_size file gel2.cluster  mode vector
fix 		chk all ave/time 10 1 4000 c_myChunk1 file T$i.density mode vector

输出文件.cluster，.density分别是液滴尺寸和密度分布了解其统计原理、输出文件的结构，得到1）最大cluster尺寸随时间变化的图 2）密度分布图。

ref: 我们自己啊

  • 用density文件得到相图。 从参考文献里可以看出，你需要根据每次温度结束时候的构象，对它进行切片处理，得到密度分布统计图，计算得到高低相的密度，用于做相图。如何切片统计密度，就需要写程序了。

    1. 读取density文件（每一帧有多少行，每一行代表什么）
    2. 得到沿着z方向的每个切片的密度图（粒子数目），类似参考文献fig 3D. 做一个差分，判断高密度相和低密度相的边界。
    3. 计算高低密度相密度。这里单位不重要。
    4. 不同帧的平均
    5. Ising拟合相图

  • 完成后回答问题：

    • 为什么要得到相图？它与实验结果会有哪些差别？
    • slab方法有哪几个优势？使用条件是什么？拉长和压缩为什么要慢？
    • density统计时候需要判断每个原子判断所处位置切片（chunk)，切片需要至少多少个？
    • 什么是系综平均？这里我们用什么来代替得到系综平均？
    • 步长选大选小有什么影响？
Statphys > Programming
#MDsimulation
Tutorial | LAMMPS 1.0
http://home.ustc.edu.cn/~lingerli/2023/01/31/lammpstutorial/
Author
Carolinge Li
Posted on
January 31, 2023
Licensed under