6.20. 复合方法 ¶


复合方法是一种复杂的脚本语言,可以直接在 ORCA 的输入中使用。使用“复合”,用户可以将正常 ORCA 计算的各个部分组合起来,以评估自定义函数。为了说明其用法,我们将使用一个示例。有关该模块的更详细描述,请参阅复合方法部分。

6.20.1. 示例


作为一个典型的例子,我们将使用描述 NH3 的“伞效应”的约束优化。该脚本将执行一系列计算,最后将打印运动的势能,并识别最小值和最大值。相应的复合脚本如下所示。

# ----------------------------------------------
# Umbrella coordinate mapping for NH3
# Author: Frank Neese
# ----------------------------------------------
variable JobName = "NH3-umbrella";
variable amin    = 50.0;
variable amax    = 130.0;
variable nsteps  = 21;
Variable energies[21];

Variable angle;
Variable JobStep;
Variable JobStep_m;
variable step;

Variable method = "BP86";
Variable basis  = "def2-SVP def2/J";

step  = 1.0*(amax-amin)/(nsteps-1);

# Loop over the number of steps
# ----------------------------
for iang from 0 to nsteps-1 do
  angle    = amin + iang*step;
  JobStep  = iang+1;
  JobStep_m= JobStep-1;
  if (iang>0) then
    Read_Geom(JobStep_m);
    New_step
      ! &{method} &{basis} TightSCF Opt
      %base "&{JobName}.step&{JobStep}"
      %geom constraints
        {A 1 0 2 &{angle} C}
        {A 1 0 3 &{angle} C}
        {A 1 0 4 &{angle} C}
        end
      end

    Step_End
  else
    New_step
      ! &{method} &{basis} TightSCF Opt
      %base "&{JobName}.step&{JobStep}"
      %geom constraints
        {A 1 0 2 &{angle} C}
        {A 1 0 3 &{angle} C}
        {A 1 0 4 &{angle} C}
        end
      end

      * int 0 1
      N 0 0 0 0.0 0.0 0.0
      DA 1 0 0 2.0 0.0 0.0
      H 1 2 0 1.06 &{angle} 0.0
      H 1 2 3 1.06 &{angle} 120.0
      H 1 2 3 1.06 &{angle} 240.0
      *
    Step_End
  endif
   Read energies[iang] = SCF_ENERGY[jobStep];
   print(" index: %3d Angle %6.2lf Energy: %16.12lf Eh\n", iang, angle, energies[iang]);
EndFor

# Print a summary at the end of the calculation
# ---------------------------------------------
print("////////////////////////////////////////////////////////\n");
print("// POTENTIAL ENERGY RESULT\n");
print("////////////////////////////////////////////////////////\n");
variable minimum,maximum;
variable Em,E0,Ep;
variable i0,im,ip;
for iang from 0 to nsteps-1 do
  angle   = amin + 1.0*iang*step;
  JobStep = iang+1;
  minimum = 0;
  maximum = 0;
  i0      = iang;
  im      = iang-1;
  ip      = iang+1;
  E0      = energies[i0];
  Em      = E0;
  Ep      = E0;
  if (iang>0 and iang<nsteps-1) then
    Em = energies[im];
    Ep = energies[ip];
  endif
  if (E0<Em and E0<Ep) then minimum=1; endif
  if (E0>Em and E0>Ep) then maximum=1; endif
  if (minimum = 1 ) then
    print(" %3d  %6.2lf %16.12lf (-)\n",JobStep,angle, E0 );
  endif
  if (maximum = 1 ) then
    print(" %3d  %6.2lf %16.12lf (+)\n",JobStep,angle, E0 );
  endif
  if (minimum=0 and maximum=0) then
    print(" %3d  %6.2lf %16.12lf    \n",JobStep,angle, E0 );
  endif
endfor
print("////////////////////////////////////////////////////////\n");

End # end of compound block


让我们先看看如何执行这个输入。为了运行它,最简单的方法是将其保存在一个普通文本文件中,使用名称“umbrella.cmp”,然后使用以下 ORCA 输入文件:

%Compound "umbrella.cmp"


不需要更多。ORCA 将读取复合文件并采取适当的行动。


关于这个 ORCA 输入的几点说明。首先,没有简单的输入行(以“!”开头)。使用复合功能时不需要简单输入,但如果用户添加了简单输入,所有来自简单输入的信息将传递给实际的复合作业。


此外,如果不想创建单独的复合文本文件,在 ORCA 中完全可以像使用其他 ORCA 块一样使用复合功能。这意味着在%Compound 指令之后,可以附加复合文件的内容,而不是给出文件名。


正如我们将看到的,在复合脚本文件中,每个复合作业可以包含正常 ORCA 输入文件的所有信息。这里有两个非常重要的例外:处理器的数量和 MaxCore。这些信息应在初始 ORCA 输入文件中设置,而不是在实际的复合文件中。


复合块的结构与所有 ORCA 块相同。它以“%”开头,以“End”结尾,如果输入不是从文件中读取的。如果复合指令在文件中,如上面的例子所示,则只需在括号内提供文件名,而不需要最后的 END。


6.20.2. 定义变量 ¶


正如我们已经指出的,可以在复合块内提供所有计算和数据处理的信息,或者创建一个包含所有细节的普通文本文件,让 ORCA 读取。后者的优点是可以使用同一个文件处理多个几何体。在前面的例子中,我们将 ORCA 指向一个外部文件。该文件“umbrella.cmp”包含所有必要的信息。


现在让我们尝试分析复合“umbrella.cmp”文件。

# ----------------------------------------------
# Umbrella coordinate mapping for NH3
# Author: Frank Neese
# ----------------------------------------------
variable JobName = "NH3-umbrella";
variable amin    = 50.0;
variable amax    = 130.0;
variable nsteps  = 21;
Variable energies[21];

Variable angle;
Variable JobStep;
Variable JobStep_m;
variable step;

Variable method = "BP86";
Variable basis  = "def2-SVP def2/J";

step  = 1.0*(amax-amin)/(nsteps-1);


第一部分包含一些一般性评论和变量定义。对于评论,我们使用与普通 ORCA 输入相同的语法,通过“#”符号。请注意,同一行中多个“#”符号会导致错误。


在初始评论之后,我们看到了一些声明和定义。声明变量的方式有很多种,详细描述见“变量 - 一般”部分。


所有变量声明以指令 Variable 开头,这表示程序将期待一个或多个新变量的声明。然后有许多选项,包括定义多个变量、为变量赋值或使用值列表。然而,所有声明必须以;符号结束。这个符号是程序的一个消息,表示当前命令的结束。在 Compound 中,每个命令末尾需要;符号是一个普遍要求,只有极少数例外。


6.20.3. 运行计算 ¶

# Loop over the number of steps
# ----------------------------
for iang from 0 to nsteps-1 do
  angle    = amin + iang*step;
  JobStep  = iang+1;
  JobStep_m= JobStep-1;
  if (iang>0) then
    Read_Geom(JobStep_m);
    New_step
      ! &{method} &{basis} TightSCF Opt
      %base "&{JobName}.step&{JobStep}"
      %geom
        constraints
          {A 1 0 2 &{angle} C}
          {A 1 0 3 &{angle} C}
          {A 1 0 4 &{angle} C}
        end
      end
    Step_End
  else
    New_step
      ! &{method} &{basis} TightSCF Opt
      %base "&{JobName}.step&{JobStep}"  
      %geom
        constraints
          {A 1 0 2 &{angle} C}
          {A 1 0 3 &{angle} C}
          {A 1 0 4 &{angle} C}
        end
      end
      * int 0 1
        N 0 0 0 0.0 0.0 0.0
        DA 1 0 0 2.0 0.0 0.0
        H 1 2 0 1.06 &{angle} 0.0
        H 1 2 3 1.06 &{angle} 120.0
        H 1 2 3 1.06 &{angle} 240.0
      *
    Step_End
  endif
  Read energies[iang] = SCF_ENERGY[jobStep];
  print(" index: %3d Angle %6.2lf Energy: %16.12lf Eh\textbackslash{}n", iang, angle, energies[iang]);
EndFor


然后我们进入信息最密集的部分。我们从 for 循环的定义开始。复合 for 循环的语法是:


对于变量 From startValue 到 endValue 执行

  指令
EndFor


如上例所示,startValue 和 endValue 可以是常数或之前定义的变量,甚至是这些变量的函数。请记住,它们必须是整数。循环结束的信号是 EndFor 指令。有关 for 循环的更多细节,请参阅 For 部分。


然后我们继续分配一些变量。

angle     = amin + iang*step;
JobStep   = iang+1;
JobStep_m = JobStep-1;


变量赋值的语法与每种编程语言中的变量相似,后面跟着 = 符号,然后是值或方程。请记住,赋值必须始终以 ; 符号结束。


下一步是每种编程语言中另一个重要部分,即 if 块。if 块的语法如下:


如果 (要评估的表达式) 那么

  指令

否则如果 ( 要评估的表达式 ) 那么

  指令
否则
  指令
EndIf


else if 和 else 部分是可选的,但最终的 EndIf 必须始终标志着 if 块的结束。有关 if 块使用的更多细节,请参阅手册的 If 部分。


接下来我们有一个特定于复合方法的命令,而不是普通编程语言的一部分。这是 ReadGeom 命令。它的语法是:


Read_Geom(整数值);


在解释这个命令之前,我们将继续处理复合脚本中的下一个命令,然后再返回到这个命令。


下一个命令是所有复合脚本的基础。这是 New_Step 命令。该命令向复合指示接下来是一个正常的 ORCA 计算。它的语法是:


New_Command 正常 ORCA 输入 Step_End


关于 New_Step 命令的一些评论。首先,在 New_Step - Step_End 命令内部,可以添加所有正常 ORCA 输入接受的可能命令。我们应该记住,定义处理器数量的命令和 MaxCore 命令将被忽略。


另一个需要记住的要点是步骤的概念。每个 New_Step - Step_End 结构对应一个步骤,从 1 开始计数(第一个 ORCA 计算)。这有助于我们定义该计算将创建的属性文件,以便我们可以用它来检索信息。


在 New_Step - Step_End 块中的一个重要特性是使用结构 &{variable}。该结构允许用户在 New_Step - Step_End 块内部使用在外部定义的变量,从而使 ORCA 输入更加通用。例如,在上面的脚本中,我们构建了计算的基本名称。

  %base "&{JobName}.step&{JobStep}"


使用定义的变量 JobName 和 JobStep。有关 &{} 结构的使用的更多细节,请参阅 & 部分,而有关 New_Step - Step_End 结构,请参阅 New_Step 部分。


最后,关于计算的几何形状有几点评论。提供几何形状给 New_Step - Step_End 计算有三种方式。第一种是传统的 ORCA 输入方式,我们可以像在所有 ORCA 输入中那样提供坐标或包含坐标的文件名。然而在 Compound 中,如果我们没有传递任何关于计算几何形状的信息,那么 Compound 将自动尝试读取上一个计算的几何形状。这是提供几何形状给复合步骤的第二种(隐式)方式。然后还有第三种方式,这就是我们在上面的例子中使用的。这是 Read_Geom 命令。该命令的语法是:

Read_Geom (步骤编号);

当我们想要将特定的几何体传递给一个计算时,可以使用此命令,该几何体并未在 New_Step - Step_End 结构中明确给出,也不是来自前一步。然后,我们只需在运行新计算之前传递我们感兴趣的计算步骤的编号。有关 Read_Geom 命令的更多详细信息,请参阅 Read_Geom 部分。


6.20.4. 数据处理 ¶


Compound 的一个最强大的功能是它可以直接访问计算的属性。为了使用这些属性,我们定义了 Read 命令。在前面的例子中,我们用它来读取计算的 SCF 能量:

Read energies[iang] = SCF\_ENERGY[jobStep];


该命令的语法是:


读取变量名 = property


其中变量名称是已经定义的变量的名称,属性是来自已知属性表中的已知属性,步骤是我们感兴趣的计算步骤。有关 Read 命令的更多详细信息,请参阅 Read 部分。

# Print a summary at the end of the calculation
  # ---------------------------------------------
  print("////////////////////////////////////////////////////////\\n");
  print("// POTENTIAL ENERGY RESULT\\n");
  print("////////////////////////////////////////////////////////\\n");
  variable minimum,maximum;
  variable Em,E0,Ep;
  variable i0,im,ip;
  for iang from 0 to nsteps-1 do
    angle   = amin + 1.0*iang*step;
    JobStep = iang+1;
    minimum = 0;
    maximum = 0;
    i0      = iang;
    im      = iang-1;
    ip      = iang+1;
    E0      = energies[i0];
    Em      = E0;
    Ep      = E0;
    if (iang>0 and iang<nsteps-1) then
      Em = energies[im];
      Ep = energies[ip];
    endif
    if (E0<Em and E0<Ep) then minimum=1; endif
    if (E0>Em and E0>Ep) then maximum=1; endif
    if (minimum = 1 ) then
      print(" %3d  %6.2lf %16.12lf (-)\textbackslash{}n",JobStep,angle, E0 );
    endif
    if (maximum = 1 ) then
      print(" %3d  %6.2lf %16.12lf (+)\textbackslash{}n",JobStep,angle, E0 );
    endif
    if (minimum=0 and maximum=0) then
      print(" %3d  %6.2lf %16.12lf    \textbackslash{}n",JobStep,angle, E0 );
    endif
  endfor
  print("////////////////////////////////////////////////////////\\n");


一旦所有数据可用,我们可以像在任何编程语言中一样在方程中使用它们。


print 语句的语法是:

print( format string, [variables]);


例如在之前的脚本中,我们这样使用它:

print(" %3d  %6.2lf %16.12lf  \n",JobStep,angle, E0 );


其中 %3d、%6.2lf 和 %16.2lf 是格式标识符,JobStep、angle 和 E0 是之前定义的变量。语法与编程语言 C 中广泛接受的 printf 命令的语法非常相似。有关 print 语句的更多详细信息,请参阅章节:Print。


与 print 命令类似的是 write2file 和 write2string 命令,它们用于写入而不是输出文件,可以选择一个文件或生成一个新字符串。


最后,确实重要的是不要忘记每个复合文件都应该以一个最终的 End 结束。


一旦我们运行前面的示例,我们将得到以下输出:

////////////////////////////////////////////////////////
// POTENTIAL ENERGY RESULT
////////////////////////////////////////////////////////
1   50.00 -56.486626696200    
2   54.00 -56.498074637200    
3   58.00 -56.505200120800    
4   62.00 -56.508823168800    
5   66.00 -56.509732863600 (-)
6   70.00 -56.508724734300    
7   74.00 -56.506590613800    
8   78.00 -56.504070086000    
9   82.00 -56.501791816800    
10   86.00 -56.500229017900    
11   90.00 -56.499674856600 (+)
12   94.00 -56.500229018100    
13   98.00 -56.501791817200    
14  102.00 -56.504070082800    
15  106.00 -56.506590613300    
16  110.00 -56.508724733100    
17  114.00 -56.509732863700 (-)
18  118.00 -56.508823172900    
19  122.00 -56.505200132200    
20  126.00 -56.498074642900    
21  130.00 -56.486626729200    
////////////////////////////////////////////////////////


随着步骤,对应步骤的角度,约束优化能量的能量加上两个极小值和一个极大值在势能中的符号。