Fork me on GitHub
image frame

蜡笔的blog

蜡笔的blog

非正交晶胞正交化方法

最近在处理的结构的时候，本来是正交的晶胞，由于我切了某个面以后，变成了非正交的了（比如说a和b并不垂直，但是c还是垂直于ab构成的平面）。懒人做法当然是上网搜教程啦~~~~，这里的内容我是从b站学习来的，有兴趣可以去原地址看视频：BV16a4y1R7qB，在B站搜索即可，我也是学习原视频学来的。

其实方法非常简单。比如说POSCAR文件，我们仅仅讨论两个不垂直的晶格矢量a和b:

这时需要考虑的是，是否存在一个线性变换，使得这两个矢量a和b可以正交。

正交的条件非常好定，只需要两个矢量a和b的余弦值为0即可（在计算机中，可以设置一个非常小的值代替0），我们现在得到条件1：

再来看看条件2。对于晶体而言，最重要的就是周期性条件，再结合一下对于POSCAR中3个矢量的理解。不难得到，变换矩阵的4个值都必须为整数，我们得到条件2：

接下来就是另一个问题，我们建系是在右手系下建立的，那么在变换以后也必须满足右手系。由行列式公式：

我们可以知道，要想保证变换后的两个矢量和c 满足右手系，这只需要这个变换矩阵的行列式为正的即可，我们得到条件3：

好了，我们已经得到这个矩阵的所有约束条件了，但是还缺一个最小化的条件4：

这样一来，我们就可以通过这样一个非线性规划问题（我没学过什么算法啊，这算是非线性规划还是线性规划啊。。。。。？？？？？？），得出来我们想要的这个矩阵。对于计算机求解，看似有点困难。但是非常幸运的是，这个矩阵中的4个数字都是整数，这样一来我们可以使用嵌套循环的方式就可以轻松解决这个问题啦！下面附Fortran代码：

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
module fun
implicit none
contains
real(kind=8) function cosin(a,b)
        real,dimension(3) :: a,b
        real :: temp1,temp2
        temp1=sqrt(a(1)**2+a(2)**2+a(3)**2)
        temp2=sqrt(b(1)**2+b(2)**2+b(3)**2)
        cosin=(a(1)*b(1)+a(2)*b(2)+a(3)*b(3))/(temp1*temp2)
end function cosin

real(kind=8) function cross(a,b)
        real,dimension(3) :: a,b
        cross=a(1)*b(2)-a(2)*b(1)
end function cross
end module fun
program nonlinear
        use fun
implicit none
integer :: x1,x2,x3,x4,s1,s2,s3,s4
integer :: i,j,k,l,m,n
integer :: up,down,temp,min_val
real,dimension(3) :: a,b,a_prime,b_prime
real :: error
a(1)=0.0  !这里从左到右输入POSCAR中不正交向量a的第1列的值
a(2)=0.0  !这里从左到右输入POSCAR中不正交向量a的第2列的值
a(3)=0.0  !这里从左到右输入POSCAR中不正交向量a的第3列的值
b(1)=0.0  !这里从左到右输入POSCAR中不正交向量b的第1列的值
b(2)=0.0  !这里从左到右输入POSCAR中不正交向量b的第2列的值
b(3)=0.0  !这里从左到右输入POSCAR中不正交向量b的第3列的值

up=10  !变换矩阵中每个整数的搜索上限
down=-10  !变换矩阵中每个整数搜索下限
error=1e-4  !最终最小化允许的误差
min_val=1000000
loop1:do i=down,up
loop2:  do j=down,up
loop3:          do k=down,up
loop4:                  do l=down,up
                                if (i==0 .and. j==0) then
                                        cycle loop4
                                end if
                                if (k==0 .and. l==0) then
                                        cycle loop4
                                end if
                                x1=i
                                x2=j
                                x3=k
                                x4=l
                                if (x2*x3-x1*x4<0) then
                                        a_prime(1)=x1*a(1)+x2*b(1)
                                        a_prime(2)=x1*a(2)+x2*b(2)
                                        a_prime(3)=x1*a(3)+x2*b(3)
                                        b_prime(1)=x3*a(1)+x4*b(1)
                                        b_prime(2)=x3*a(2)+x4*b(2)
                                        b_prime(3)=x3*a(3)+x4*b(3)
                                        !write(*,*) abs(cosin(a_prime,b_prime))
                                        if (abs(cosin(a_prime,b_prime))<=error) then
                                                temp=abs(x1)+abs(x2)+abs(x3)+abs(x4)
                                                if (temp<min_val) then
                                                        min_val=temp

                                                        s1=x1
                                                        s2=x2
                                                        s3=x3
                                                        s4=x4
                                                end if
                                        end if
                                end if
                        end do loop4
                end do loop3
        end do loop2
end do loop1
write(*,*) s1,s2,s3,s4,min_val
stop
end program

别看套了4个循环，效率还可以。需要自己改误差参数以及矢量的分量值。输出的结果右4个数字，就是x1,x2,x3,x4的值，最后一个值是最小化的值，用来判断是否出错。

改完以后就可以用M某软件重新定义晶胞操作（具体可以Google一下怎么做），即可得到正交后的晶胞（但是不是立方，仅仅只能正交）。对于三个矢量都不垂直的情况，可以再探究探究。需要注意的是，重新定义后的正交晶胞的体积会变大，M某软件会提示你体积变大了多少倍，所以最后还要自己看情况调整。

最后最重要的是，任何变幻出来的正交胞，需要自己再三确认是否是原来的结构！！！！！！！！！！！！！！这非常重要，可以使用XRD辅助检验，一定要手动检查。

懒人最近的施工计划：

  1. KNN最近邻算法

  2. 基于高斯核函数的密度估计算法

  3. 晶体空间群识别算法（这个论文还没看完捏~~）

就这样啦，拜拜~~~~~~~~~~~~~~~~~~~~

最近写的一些全自动处理脚本

最近写了一个扩胞全自动的脚本，使用Fortran编写，挂上来以防哪天丢了。

第1个是一个扩胞的脚本，支持的是POSCAR格式。下载链接：GitHub - whitecrn/crystal_enlarge文件中的原子支持前111号原子自动识别（包括第111号原子）；7种晶系全都支持；分数坐标和绝对坐标会自动识别。单个结构中除了不要出现112号以后的元素以外，没什么大的限制。但是，文件中不能出现原子电荷信息！！！！比如Li+、F-这种。

使用方法：需要有gnu库的支持。将文件与需要扩胞的文件放在同一目录下。使用：

1
chmod +x crystal_enlarge

然后再使用命令：

1
./crystal_enlarge (x方向扩胞倍数) (y方向扩胞倍数) (z方向扩胞倍数) (需要扩胞的POSCAR文件名) (扩胞后输出的文件名)

举个例子：

1
./crystal_enlarge 2 3 4 POSCAR POSCAR-ENLARGE

上述代码就是将POSCAR文件进行扩胞操作，在x方向上×2 ， y方向上×3，z方向上×4，并且将扩胞的结果输出到POSCAR-ENLARGE里面。

经过本人测试，从material project上随意下了几个三斜和单斜的晶胞，扩胞出来的结果都没有问题，略略略略略略！！！有技术问题或者要源代码可以想办法联系我~~0.0~

哦对了，如果有需要转换以后的POSCAR格式文件转换成lammps文件的格式。可以参考我之前的一个poscar2data的文件脚本。所支持的体系也和上面的类似，不支持电荷信息，支持前111号原子的自动识别（包括111）转换。参考链接：VASP的POSCAR格式转换lammps格式脚本 | 蜡笔的blog。关于库的支持，在poscar2data这个网页里面有很详细的讲解。

下面举一些例子（随便在material project上找的结构）。和vesta软件扩出来的完全没有区别，但可以进行超大规模扩胞（windows下的vesta扩太大直接会变得很卡），以及命令行批量操作。 不过批量使用之前，最好先用vesta扩一个胞，经过对比无误后再进行大规模扩胞（手动保命）。

我扩的三斜（mp-1176991, Li6MnSb3(PO4)6）：

vesta扩的三斜（mp-1176991）：

我扩的单斜（mp-705035, Mo(PO3)4）:

vesta扩的单斜：

我扩的六方（mp-1227397, BiPt2Pb）：

vesta扩的六方：

我扩的三方（mp-1227397, CsGeCl3）：

vesta扩的三方：

我扩的正交（mp-1227397, LiFe7O7F）：

vesta扩的正交：

Fortran快速排序

处理数据的时候经常要用排序法，但是无可奈何我好像找不到gfortran上现成的排序程序。于是自己在写数学函数库的时候，也加入了这一部分的内容。

起先写的是冒泡排序，但是同学建议我冒泡排序太慢了，对于数量多起来的排序时间会变得非常夸张，于是乎3.29下午开始学习快速排序。3.30中午终于把gfortran的快速排序代码写好了。。。。。 算法详解什么的有空再写吧，先贴代码。。。。

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
recursive subroutine Fast_Sort(A,sort_type)
        real :: A(:)
        character(len=*),optional :: sort_type
        integer :: head,tail,i,s
        real :: temp
        real :: pivot
        if (.not. present(sort_type)) then
                sort_type='asc'
        end if
        head=0
        tail=size(A)+1
        pivot=A(1)
loop1:  do while (.true.)
                tail=tail-1
                head=head+1
tail_loop:      do while (.true.)
                        if (A(tail) <= pivot) then
                                exit tail_loop
                        end if
                        tail=tail-1
                end do tail_loop

head_loop:      do while (.true.) 
                        if (A(head) > pivot) then
                                exit head_loop
                        end if
                        head=head+1
                        if (head > size(A)) then
                                exit head_loop
                        end if
                end do head_loop

                if (head < tail) then
                        temp=A(head)
                        A(head)=A(tail)
                        A(tail)=temp
                else if (head == tail) then
                        temp=A(head)
                        A(head)=A(1)
                        A(1)=temp
                        exit loop1
                else if (head > tail) then
                        temp=A(tail)
                        A(tail)=A(1)
                        A(1)=temp
                        exit loop1
                end if
        end do loop1
                if (size(A(:tail-1)) > 1) then
                        call Fast_Sort(A(:tail-1))
                end if

                if (size(A(tail+1:)) > 1) then
                        call Fast_Sort(A(tail+1:))
                else 
                        return
                end if

        if (trim(adjustl(sort_type)) == 'desc') then
                call reverse_array(A)
        end if


contains

subroutine reverse_array(array)
        implicit none
        real :: array(:)
        integer :: s,i,temp

        s=size(array)

        do i=1,s/2
                temp=array(i)
                array(i)=array(s+1-i)
                array(s+1-i)=temp
        end do
end subroutine reverse_array

end subroutine Fast_Sort

使用了Fortran递归的方法写的快速排序算法。其实本体只有升序排序，如果可选输入参数为‘desc’，这个时候就会启动后面的reverse_array程序，将升序排号的结果翻转就是降序啦（偷懒第一名。。。。，其实是不会写降序，不要问我为什么？？？）。子程序我是放在模块里面使用的，所以实际上需要给这个子程序增加一个显式接口，关于显式接口可以上网搜。简单暴力操作就是在这个程序上下方加上如下内容：

1
2
3
4
5
6
module xxxx(自己定模块的名字)
contains

#这里放上面的代码

end module xxxx(和前面的名字一致)

如果想要将子程序和主程序放在一起，可以将上述内容放在主程序之前(模块必须放在主程序之前)然后在主程序中，在implicit none之前写上：

1
2
use xxx(这个名字是module的名字)
implicit none

然后在需要使用的时候就可以：

1
call Fast_Sort(A,'asc/desc')

A为需要排序的数组，后一个参数可以输入字符串’asc’和’desc’，输入’asc’的时候为升序排列，’desc’时为降序排列，**当不输入第二个参数或者输入的第二个参数不是‘asc’和’desc’其中之一时，默认为升序。**

最后是对于10w个随机数，使用冒泡排序和快速排序二者所耗时间比较，时间是cpu时间，单位应该是s吧？

效果还是很明显的，快速排序耗时远远小于冒泡排序，嘻嘻~~

至于如果想把子程序和主程序分开放，可以上网搜一下Fortran的静态库或者动态库的编译和链接，这里就不详细讲了，认识我的话可以来问我hhhhh~

就这样啦，拜拜~~

VASP的POSCAR格式转换lammps格式脚本

近期处理POSCAR转换lammps的data文件的时候，需要ovito协助处理，但是命令行的大量处理好像还是成问题的（特别是我几百个POSCAR需要转换的时候，一个个手拖就不现实了）

之前写过一些自动转换格式的脚本，但是只是针对于某一个特定的结构进行的，多一个原子少一个原子，原子换一下就用不了的。而我每次处理的POSCAR文件的原子种类和个数都在发生变化。于是一怒之下就写了这个通用的脚本，可以直接将含有前111号原子的POSCAR文件自动转换成lammps可读取的data文件格式。*支持分数坐标和绝对坐标的读取。下载链接：GitHub - whitecrn/poscar2data: This scrip can help you transform POSCAR file to lammps data file

 本脚本使用gfortran编写，使用的是gnu编译器。脚本所需要的依赖如下图所示。

 文件的源代码我没有放上来，直接上传的是编译好的文件。因为源代码里面有很多我自己写的库文件，最终合成打包起来比较麻烦，有需要的话可以联系我问我拿。源代码各文件如下图。。。。。。

POSCAR需要比较固定的开头格式，中间不能有空行（如果格式不是固定的，ovito也读取不了。。。。）。格式图下图所示。

第一行是体系的名称，第二行是缩放因子，第三行到第五行是晶格参数，第六行是原子种类信息，第七行是原子对应的原子数，第八行是区分绝对坐标和分数坐标，第九行以后开始写的是原子位置信息。目前只能做到基于此类的转换，如果含有电荷信息可能无法使用。原子的种类限制在前111号元素，原子的个数和种类个数都不会受到限制。

使用脚本很简单，下载以后（Linux系统上需要装有gnu编译器）。将脚本放在与需要处理的POSCAR文件同一目录之下，使用命令：chmod +x poscar2data,为脚本增加可执行权限。

然后如果需要转换的POSCAR格式文件名称就是POSCAR，那使用命令：**./poscar2data**，将会自动转换POSCAR文件，输出成一个新的data文件：lammps.data，即可完成。

 如果你的需要转换的文件不叫POSCAR，需要输出的文件想换一个名称，可以使用命令：**./poscar2data (输入文件名称) (输出文件名称)** ，例如: ./poscar2data POSCAR1 FeCoNi.data，那么脚本就会读入POSCAR1文件，输出lampps格式的FeCoNi.data文件。

！！！注意，如果你要输入可选参数的时候，第一个必须为输入文件的名称，第二个必须为输出文件的名称，如果只输入一个文件名称，则默认这个名称为输入文件的名称。

转载请注明出处：GitHub whitecrn

一些计划和想法

发现寒假要做的事情太多了，书也看不过来了，大概列一下寒假想要学的东西的先后顺序吧。学是肯定学不完的，学完了中文的书还要看经典著作，太难了我。

  1. 《热力学与统计物理》。这本感觉我必须要在寒假学完，MD很多内容（系综之类的，蒙特卡洛模拟啥的，还有固体物理很多分布）都和热统有关系。本科的时候只自学过《热学》，在赵凯华先生的书上顺带学了一些统计物理的内容，但是不够全面。我是想优先使用汪志诚先生的《热力学·统计物理》，这本内容看起来并不是非常难，我觉得可以很快的扫过。还准备了朗道的《统计物理学 Ⅰ》，朗道这本热统感觉比其他九卷看起来难度要低一些，我觉得学完汪先生的可以再看看朗道的。

  2. 《固体物理》。这本没啥好说的，第二重要。先把胡安先生的课跟完，再复习推导一下前面几章的内容，然后有时间的话看一下ASHCROFT/MERMIN两位的《SOLID STATE PHYSICS》。

  3. 《电动力学》。这个蛮重要的，量子场论里面提到了非常多这个内容，直接跟着田光善老师的课就好了。其实我一直想看王振林老师的《现代电动力学》，但无可奈何时间不够。

  4. 《量子场论》。这个课挺有意思的，但是帮助不是非常非常大，所以我排在了第四个。打算直接啃Pekin的《An Introduction to Quantum Field Theory》，能学多少是多少吧这个。

  5. 《Fortran》。没啥好说的，希望把《Fortran 95程序设计》剩下的内容补完，这个应该不会花太长时间。

  6. 《机器学习》。这个主要是大概了解整个机器学习原理和过程，懂得分析数据就行，中间的推导和具体公式应用就不花时间学了，不然是无止境的。

  7. 《量子力学》。《量子力学》和《高等量子力学》都学完了，但是还是希望有空可以看一下科恩·塔诺季的《量子力学》。我觉得这本书看到现在为止，很多内容分析和讲的都非常透彻，感觉可以多读读。

  8. 剩下就是一些计算机上的一些计划，不过真正做完都不知道啥时候了。想把量子力学学到的公式打一个总结文档出来；编写一个Fortran子程序或者库，把这学期《数值分析》学到的各种算法写成Fortran，这样以后要使用或者调用也方便（已经有matlab代码了，翻译一下就好了）；还有想把科研过程中学到的各种数据处理之类或者经验上的内容分享到博客上来，这个看有时间就更新吧，差不多就这样。

  • Copyrights © 2015-2024 白色的蜡笔
  • Visitors: | Views:

请我喝杯咖啡吧~

支付宝
微信