0

我有水分子速度的模拟数据。数据格式如下。为了清楚起见，我想描述数据的格式，它很容易导致我想要计算的内容。

一个水分子由三个原子组成：氧（O）和两个氢（H）。在这里，我将它们命名为 O、H1 和 H2。

下面的数据以 linetitle 0和 number开头4335，表示它包含 4335 个原子（4335/3 = 1445 个水分子）。

从第三行 ( 0.0923365 0.0341984 -0.1248516) 开始的前三个数字表示氧 (O) 原子在三个笛卡尔方向 Ox、Oy、Oz 上的速度。接下来的三个数字，在同一行中代表氢 (H1) ==> H1x、H1y、H1z 的速度。最后，第四行中的前三个数字代表氢 (H2) ==> H2x,H2y,H2z 的速度。最后，同一第四行中的以下三个数字代表氧原子的速度。

这些序列对 2170 行中的所有 4335 个原子进行，包括数据文件中的前两行，并在下一节从 开始重复title 1。

title 0
4335  2.0001000e+04
 0.0923365   0.0341984  -0.1248516  -0.8946258   1.6688854   0.8259304
 0.2890579   0.8051153  -1.5612963   0.0625492  -0.1361579   0.2869132
 0.2343408  -0.0665305   1.0745378  -0.8375892   0.6953992   0.5149021
-0.1628550   0.0131844   0.0688080   0.2429340   0.2168210  -0.0289806
-0.3677613   0.2054004  -0.1511643  -0.3487551  -0.1454157   0.0801884
-0.9039297  -0.0682939  -0.2337404  -0.5605327  -0.0369157   0.2243892
-0.3100274  -0.2673132  -0.2093299   0.1975043  -0.4572202  -0.8410826
-0.6995287  -0.4123909   0.0649209  -0.1910519   0.2289656   0.2443295
-0.0279093   0.5790939  -0.0104249  -1.1961776  -0.5387340   0.1445187
-0.3188485   0.3789352  -0.0112114   0.7831523   0.6043882  -0.7131590
-0.7214440  -0.5358508  -0.3035673  -0.1549275  -0.1402387  -0.0101964
-0.2027608   1.5107149   0.2963312  -1.5104872  -0.1554981  -1.3323215
 0.1097982  -0.1553742   0.3803437   0.0816858   0.0265007   0.4215823
 0.1157368   0.2100116   0.4712551   0.1799426  -0.1260255  -0.2131755
 0.1811777  -0.9442581  -0.6036636   0.9681703  -0.1523646  -0.3502441
 0.0976771   0.0019619  -0.1832204  -0.0055989   0.2701100  -0.4416720
 0.8496723   0.4070951  -0.0819204   0.1156806  -0.1619873  -0.0016126
-0.4051959   0.4263505  -0.9460036   0.4412067   0.1002270   0.5864405
-0.3831136   0.3240860  -0.0005143  -0.5667163   0.2618876   0.0103317
-0.6442209   0.3965833  -0.0778050  -0.2404238  -0.1339887  -0.1662417
 0.3421198   0.7480828  -1.8316993  -0.4454920  -0.0499657  -0.1951254
-0.2895359  -0.1934811  -0.2674928   0.1255802   1.3522828  -0.2829485
-0.4129106  -0.6842645  -1.0147657  -0.1278501  -0.0597648  -0.1478294
-0.2519974   0.0665314  -0.0690079  -0.0480210  -0.1179547  -0.2091919
-0.1942484   0.2583650  -0.0734658  -0.1216313   0.5158040  -0.0676843
-0.3063602   0.8148463  -0.1959571  -0.1009838  -0.3394633  -0.0866587
 .
 .     (goes on until line 2170)
 .
 0.1028815  -0.0844088  -0.2156557  -0.1698745  -0.2018967  -0.3863209
 0.1793070  -0.1005802   0.1800752  -0.1404713   0.2216020   0.2236271
 0.5192825  -0.7398186   0.0418758   0.0347715  -0.3457840  -0.1300237
-0.3089482   1.1125441  -0.4020403   0.2739744  -0.9062766   0.0012294
 0.1498538   0.0883857  -0.0094638   0.0963565  -1.1027019   0.0115313
-0.0432824   0.3330713   0.0304943
title 1
4335  2.0002000e+04
-0.2082078   0.1774843  -0.1023302  -0.1100437   0.5973607   1.0627041
-0.2216015   0.0448885  -0.8415924   0.1691296   0.6008261  -0.0373434
 0.9387534  -0.3642305   0.6756270  -0.6000357   0.6632088   1.0567899
-0.3234407  -0.1781680  -0.1936070  -0.4799916  -0.1522612  -0.2347461
 0.1045985   0.1999704  -0.1482928  -0.0439331   0.0413923   0.1605458
 0.3403952  -0.2012104   0.4851457  -0.9665228   0.2202362   0.0046218
 .
 .   (goes on until line 2170)
 .

我想计算的是每个分子的合成速度，我想用 Perl 来做这个。算法是这样进行的。

首先将氧气 (O) 和氢气 (H1 & H2) 的速度分别存储在 Ox,Oy,Oz, H1x,H1y,H1z 和 H2x,H2y,H2z 中。

接下来我们定义：

    velocity_x = Ox + Hx + Hx
    velocity_y = Oy + Hy + Hy
    velocity_z = Oz + Hz + Hz

最后计算

  resultant_velocity = sqrt(velocity_x**2 + velocity_y**2  + velocity_z**2)

并将“resultant_velocity”存储到新文件中（文件应该是 title_0.dat）。并且程序将计算文件中从title 1直到开始的速度title 200。

我是 Perl 的新手，但我想在 Perl 中做这个操作，因为我发现它非常有趣。我可以在 Perl 中编写简单的“读写”操作，但不知道如何拆分并将值分配给变量并执行计算，尽管计算是高中标准。

#!/usr/bin/perl -w

$data_file="malto.dat";
open(DAT, $data_file) || die("Could not open file!");
@raw_data=<DAT>;
close(DAT);

while(<@raw_data>){
@columns=split /\s+/,$_; 
if($columns[0]=~ m/ATOM/){ 
print "$columns[5], $columns[6], $columns[7]\n";
    }
}

我想从专家那里得到一些指导，这样我可以在编写代码的同时加深对 Perl 的理解。

感谢任何帮助。问候

4

3 回答 3

1

这是我的建议：将工作分解为小组件，并为工作的每个有意义的部分编写一个方法。以机智：

use strict;
use warnings;

main(@ARGV); # Pass data file name on command line. Don't hard-code it.

sub main {
    my $data_f = shift;
    open(my $data_h, '<', $data_f) or die "$!: $data_f";
    while (my $section = get_section($data_h)){
        # Also write methods that can be called here to make 
        # desired computations, print output, etc.
    }
}

sub get_section {
    # Takes a file handle.
    # Returns a hash reference containing all of the data
    # for an entire section of the file.
    my $h = shift;
    return if eof($h);
    chomp (my $title = <$h>);
    my ($n_atoms) = <$h> =~ /^(\d+)/;
    return {
        'title'     => $title,
        'n_atoms'   => $n_atoms,
        'molecules' => get_molecules($h, $n_atoms / 3),
    };
}

sub get_molecules {
    my @molecules;
    return \@molecules;
}

我没有写get_molecules()方法。它需要一个文件句柄和一个整数（N 个分子）。它可以返回对数组数组或哈希数组的引用，每个内部数组/哈希都保存单个分子的信息。

于 2012-12-14T15:57:09.880 回答
1

也许以下内容会对您有所帮助：

use strict;
use warnings;
use Math::Complex;

my $dataFile = 'malto.dat';

{
    local $/ = 'title ';
    open my $fh, '<', $dataFile or die $!;

    while (<$fh>) {
        chomp;

        my @data = split or next;

        my $titleNum     = 'Title ' . shift @data;
        my $atom         = shift(@data) . ' ' . shift @data;
        my $resultantVel = calcResultantVel( \@data );

        print $titleNum, "\n";
        print $atom,     "\n";
        print 'ResultantVel: ' . $resultantVel, "\n\n";
    }

    close $fh;
}

sub calcResultantVel {
    my ($dataRef) = @_;

    my ($velocity_x, $velocity_y, $velocity_z);

    while ( my @nums = splice( @$dataRef, 0, 9 ) ) {
        $velocity_x += $nums[0] + $nums[3] + $nums[6];
        $velocity_y += $nums[1] + $nums[4] + $nums[7];
        $velocity_z += $nums[2] + $nums[5] + $nums[8];
    }

    return sqrt( $velocity_x**2 + $velocity_y**2 + $velocity_z**2 );
}

单词和空格title用作记录分隔符，因此每次读取都会接收一块由 . 分隔的数据title。chomp删除记录分隔符，然后记录在split空白处。

第零个元素是标题号，这是shifted off @data。的第一个和第二个元素@data是原子数，它们也被shift删除了。剩余的数组元素是笛卡尔方向，并且对该数组的引用被发送到子例程calcResultantVel。

该子例程一次采用九个元素的块：三个用于 O 原子，三个用于第一个 H 原子，三个用于第二个 H 原子，并且根据您提供的定义保留运行总和。最后，返回结果速度。

这是一些示例输出：

Title 0
4335 2.0001000e+04
ResultantVel: 13.2945751170603

Title 1
4335 2.0001000e+04
ResultantVel: 12.7696611061461

您可以直观地验证它是否正常工作。由于您“……可以在 Perl 中编写简单的‘读写’操作……”，下一步就是让它将所需的结果写入文件。

希望这可以帮助！

于 2012-12-14T22:12:51.783 回答
0

感谢您的帮助和指导。我已尝试修改您的代码，如下所示。它至少可以满足我的需要。

#!/usr/bin/perl
###############

#use strict;
#use warnings;
use Math::Complex;

open    OUTPUT, '>', "velocityOnly.dat" or die "Can't create filehandle: $!";

my $dataFile = 'velF1F2.vel';

{
local $/ = 'title ';
open my $FH, '<', $dataFile or die $!;

while (<$FH>) {
    chomp;

    my @data = split or next;

    my $titleNum     = 'Title ' . shift(@data);
    my $atom         = shift(@data) . ' ' . shift(@data);
    #my $resultantVel = calcResultantVel( \@data );

    #print OUTPUT "$titleNum", "\n";
print "$titleNum", "\n";

for my $i (1..1445)
{
    $j=(9*($i-1));  

   $velocity_x = $data[($j+0)] + $data[($j+3)] + $data[($j+6)];
   $velocity_y = $data[($j+1)] + $data[($j+4)] + $data[($j+7)];
   $velocity_z = $data[($j+2)] + $data[($j+5)] + $data[($j+8)];

   $velo = sprintf '%.3f',sqrt( $velocity_x**2 + $velocity_y**2 + $velocity_z**2 );
        chomp $velo;
        print   "$velo","\n";
        print   OUTPUT  "$velo\n";

    }

    #print 'ResultantVel: ' . $resultantVel, "\n\n";
    }

   close $FH;
   }

但我想通过添加一些其他功能来进行一些复杂的计算来进一步扩展。编码

在此之前，我需要一些关于将以下代码变成子程序的指南。我在这里有点迷路了。您的代码实际上添加了所有 X、Y 和 Z，最后找到了速度。但我想要的不是那个。每个 9 个值随后代表包含三个原子的水分子的坐标。

（数字 1445 是分子数。每个分子包含三个原子，每个原子有三个坐标。所以对于一个水分子有 9 个笛卡尔坐标。）

这里的 i 代表水分子的数量 

for my $i (1..1445)
{
    $j=(9*($i-1));  

   $velocity_x = $data[($j+0)] + $data[($j+3)] + $data[($j+6)];
   $velocity_y = $data[($j+1)] + $data[($j+4)] + $data[($j+7)];
   $velocity_z = $data[($j+2)] + $data[($j+5)] + $data[($j+8)];

   $velo = sprintf '%.3f',sqrt( $velocity_x**2 + $velocity_y**2 + $velocity_z**2 );
        chomp $velo;
        print   "$velo","\n";
        print   OUTPUT  "$velo\n";

    }
于 2012-12-17T12:45:06.950 回答