Chain Alignment Parameter （链取向参数）

关注 $Mr.material,$ $更$$多$$精$$彩$！

---

主要专栏内容包括：
$\text{ †《LAMMPS小技巧》：}$ 主要介绍采用分子动力学（$Lammps$）模拟相关安装教程、原理以及模拟小技巧（难度：$★$）
$\text{ ††《LAMMPS实例教程—In文件详解》：}$ 主要介绍采用分子动力学（$Lammps$）模拟相关物理过程模拟。（包含：热导率计算、定压比热容计算，难度：$★$$★$$★$）
$\text{ †††《Lammps编程技巧及后处理程序技巧》：}$ 主要介绍针对分子模拟的动力学过程（轨迹文件）进行后相关的处理分析（需要一定编程能力。难度：$★$$★$$★$$★$$★$）。
$\text{ ††††《分子动力学后处理集成函数—Matlab》：}$ 主要介绍针对后处理过程中指定函数，进行包装，方便使用者直接调用（需要一定编程能力，难度：$★$$★$$★$$★$）。
$\text{ †††††《SCI论文绘图—Python绘图常用模板及技巧》：}$ 主要介绍针对处理后的数据可视化，并提供对应的绘图模板（需要一定编程能力，难度：$★$$★$$★$$★$）。
$\text{ ††††††《分子模拟—Ovito渲染案例教程》：}$ 主要采用 $\mathrm{Ovito}$软件，对 $\mathrm{Lammps}$ 生成的轨迹文件进行渲染（难度：$★$$★$）。

$\text{ 专栏说明（订阅后可浏览对应专栏全部博文）：}$
$注意：$如需只订阅某个单独博文，请联系博主邮箱咨询。$lammps\_materials@163.com$

$♠$$†$ 开源后处理集成程序：请关注专栏《LAMMPS后处理——MATLAB子函数合集整理》
$♠$$†$$†$ 需要付费定制后处理程序请邮件联系：$lammps\_materials@163.com$

---

转载一个Chain Alignment Parameter （链取向参数）的计算code

#This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.


### Purpose: This Script will compute Hermans orientation factor (between e2e vector and direction of pull) for systems undergoing deformation ####
### Syntax: python orientation.py < (filename).dump > orientation.txt ###
### Author: Janani Sampath ###
### Date: Feb 2016 ###


#This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.


import sys, string
from numpy import *
from math import *
import numpy as np
import math as math

tframes = 1
skiptoframe = 0
linelist = []
coslist = [0]*300
natoms = 360972
nbeads = 160
npoly = natoms/nbeads
strain = [0]*(300)
f = [0]*(300)
count = [0]*300

for i in range(0, skiptoframe):
    sys.stdin.readline() 
    sys.stdin.readline()
    sys.stdin.readline()
    sys.stdin.readline()
    sys.stdin.readline()
    sys.stdin.readline()
    sys.stdin.readline()
    sys.stdin.readline()
    sys.stdin.readline()

    for j in range(natoms):
        line = sys.stdin.readline()      
  
for i in range(0, tframes):
    sys.stdin.readline() 
    sys.stdin.readline()
    sys.stdin.readline()
    sys.stdin.readline()
    sys.stdin.readline()
    line = sys.stdin.readline()
    [xlow,xhi] = map(float,line.split())
    xbox = xhi-xlow
    line = sys.stdin.readline()
    [ylow,yhi] = map(float,line.split())
    ybox = yhi-ylow
    line = sys.stdin.readline()
    [zlow,zhi] = map(float,line.split())
    zbox = zhi-zlow
    sys.stdin.readline()

    for j in range(natoms):
        line = sys.stdin.readline()
        [ii,molj,typej,xu,yu,zu] = string.split(line)
        typemon = int(typej)
        w=int(ii)
        x = float(xu)  ##unscale and unwrap coordinates
        y = float(yu)
        z = float(zu) 
        line=[w,x,y,z,xbox,typemon]
        linelist.append(line)
           
for i in range(0,natoms-1):
    if linelist[i][5] == 3 or linelist[i+1][5] == 1:
        continue
    else:
        x1 = linelist[i][1]
        x2 = linelist[i+1][1]
        y1 = linelist[i][2]
        y2 = linelist[i+1][2]
        z1 = linelist[i][3]
        z2 = linelist[i+1][3]
        ABy = (y2-y1)           #bond vector coordinates 
        if abs(ABy) > 5:        #check if two beads are straddling the periodic boundary (check distance must be some number larger than 1.5 which is maximum extent of bond)
            continue            #do not include these bond vectors in the calculation
        else:
            avgloc = int((y2+y1)/2)
            count[avgloc] += 1
            ABx = (x2-x1)               #bond vector coordinates     
            if abs(ABx) > 5:            #check if two beads are straddling the periodic boundary
                if ABx < 0:
                    ABx = ABx + xbox 
                elif ABx > 0:
                    ABx = ABx - xbox
            ABz = (z2-z1)                              #bond vector coordinates
            AB = math.sqrt((ABx**2)+(ABy**2)+(ABz**2)) #bond vector between two beads
            ABdotBC = ABz*zbox                         #dot product of bond vector with (x, y, or z box coordinate)
            cosT = (ABdotBC)/((AB)*zbox)               #cosine theta (angle between bond vector and box coordinate)
            cosT2 = cosT*cosT                          #squared
            coslist[avgloc] += cosT2                   #sum of all values in (avgloc) bin

for i in range(0, len(coslist)):
    if count[i] == 0:
        continue
    else:
        mean = coslist[i]/count[i]
        f[i] = ((3*mean)-1)/2

print "strain orientation"
for k in range(0,len(f)):
    print k,f[k]