ORCA 예제는 ORCA 가이드에 실린 ORCA 6 출력이다. opt.out은
물의 B3LYP-D4/def2-TZVP 구조 최적화 + 진동수 계산이다. QE와 LAMMPS 예제도 각 가이드를 쓰면서 실제로 돌린 출력이다.
grep, sed, awk 문법은 앞 장에 정리했다.
ORCA: 마지막 에너지
ORCA는 에너지를 계산할 때마다 FINAL SINGLE POINT ENERGY 줄을 쓴다. 최적화에서는 사이클마다 한 줄씩 나오므로
마지막 줄만 가져와야 한다.
$ grep 'FINAL SINGLE POINT ENERGY' opt.out
FINAL SINGLE POINT ENERGY -76.426286364586
FINAL SINGLE POINT ENERGY -76.426353431114
FINAL SINGLE POINT ENERGY -76.426364583348
FINAL SINGLE POINT ENERGY -76.426366351564
FINAL SINGLE POINT ENERGY -76.426369457144
$ grep 'FINAL SINGLE POINT ENERGY' opt.out | tail -1 | awk '{print $5}'
-76.426369457144
awk 하나로도 된다. 맞는 줄마다 변수에 덮어쓰고 END에서 출력하면 마지막 값이 남는다.
줄이 하나도 없으면 빈 줄이 나오므로, 스크립트에서는 결과가 비었는지 확인한다.
$ awk '/FINAL SINGLE POINT ENERGY/ {e = $5} END {print e}' opt.out
-76.426369457144ORCA: 최적화 사이클별 에너지
수렴 과정을 그림으로 보려면 사이클 번호와 에너지를 두 열로 저장한다. gnuplot이나 matplotlib이 그대로 읽는다.
$ grep 'FINAL SINGLE POINT ENERGY' opt.out | awk '{printf "%d %.10f\n", NR, $5}' > scf.dat
$ cat scf.dat
1 -76.4262863646
2 -76.4263534311
3 -76.4263645833
4 -76.4263663516
5 -76.4263694571
다섯 줄은 최적화 네 사이클과 수렴한 구조에서의 마지막 계산이다. 마지막 값 기준 상대 에너지(kcal/mol)를 셋째 열로 붙이면
수렴이 얼마나 남았는지 한눈에 보인다. 모든 값을 배열에 모은 뒤 END에서 계산한다.
$ awk '/FINAL SINGLE POINT ENERGY/ {n++; e[n] = $5}
END {for (i = 1; i <= n; i++) printf "%d %.10f %8.4f\n", i, e[i], (e[i] - e[n]) * 627.5095}' opt.out
1 -76.4262863646 0.0521
2 -76.4263534311 0.0101
3 -76.4263645833 0.0031
4 -76.4263663516 0.0019
5 -76.4263694571 0.0000ORCA: 정상 종료와 수렴 확인
값을 쓰기 전에 계산이 끝까지 돌았는지부터 본다. 정상 종료한 ORCA 출력의 마지막 부분에는
ORCA TERMINATED NORMALLY가 있고, 최적화가 수렴하면 THE OPTIMIZATION HAS CONVERGED가 찍힌다.
05-killed는 비교를 위해 출력 앞부분만 남겨 walltime에 걸려 죽은 계산처럼 만든 파일이다.
$ grep -c 'ORCA TERMINATED NORMALLY' opt.out
1
$ grep -c 'ORCA TERMINATED NORMALLY' jobs/05-killed/water.out
0
$ grep -q 'THE OPTIMIZATION HAS CONVERGED' opt.out && echo converged || echo 'NOT converged'
converged
grep -q는 아무것도 출력하지 않고 종료 코드만 돌려준다. 찾으면 0, 못 찾으면 1이라 if에 바로 쓴다.
ORCA: 허수 진동수와 자유 에너지
진동수는 VIBRATIONAL FREQUENCIES와 NORMAL MODES 사이에 있다. awk의 범위 패턴
/시작/,/끝/으로 그 구간만 잘라 낸 뒤 둘째 열을 본다.
$ awk '/VIBRATIONAL FREQUENCIES/,/NORMAL MODES/' opt.out | awk '/cm\*\*-1/ {print $2}'
0.00
0.00
0.00
0.00
0.00
0.00
1615.68
3780.78
3885.42
처음 여섯 개의 0은 병진과 회전이다. 허수 진동수는 ORCA가 음수로 쓰고 뒤에 ***imaginary mode***를 붙인다.
아래 명령은 구간 안에서 음수인 줄만 세어 출력한다. ts.out은 형식을 보이려고 6번 모드를
ORCA 가이드의 예시처럼 음수로 고친 사본이다.
$ awk '/VIBRATIONAL FREQUENCIES/ {f = 1} /NORMAL MODES/ {f = 0}
f && /cm\*\*-1/ && $2 < 0 {n++; print}
END {print n + 0, "imaginary"}' opt.out
0 imaginary
$ awk '/VIBRATIONAL FREQUENCIES/ {f = 1} /NORMAL MODES/ {f = 0}
f && /cm\*\*-1/ && $2 < 0 {n++; print}
END {print n + 0, "imaginary"}' ts.out
6: -456.78 cm**-1 ***imaginary mode***
1 imaginary
n + 0은 한 번도 증가하지 않은 변수를 빈 문자열이 아닌 0으로 출력하게 한다.
열화학 결과는 이름이 붙은 줄 하나라 grep과 awk로 바로 뽑힌다.
$ grep 'Final Gibbs free energy' opt.out
Final Gibbs free energy ... -76.42286885 Eh
$ awk '/Final Gibbs free energy/ {print $6}' opt.out
-76.42286885ORCA: 마지막 구조를 xyz로
CARTESIAN COORDINATES (ANGSTROEM) 블록은 사이클마다 나온다. 머리글 다음 줄(-----)을 건너뛰고
빈 줄이 나올 때까지 좌표를 모은다. 새 블록을 만나면 개수를 0으로 되돌리므로 끝에는 마지막 블록만 남는다.
awk '/CARTESIAN COORDINATES \(ANGSTROEM\)/ {getline; n = 0; blk = 1; next}
blk && NF == 0 {blk = 0}
blk {line[++n] = $0}
END {print n; print "final geometry"; for (i = 1; i <= n; i++) print line[i]}' opt.out > final.xyz$ cat final.xyz
3
final geometry
O 0.000000 0.000000 0.111489
H 0.000000 0.765189 -0.473160
H 0.000000 -0.765189 -0.473160
정규식 안의 괄호는 \(로 이스케이프한다. getline은 다음 줄을 읽어 버리는 명령이다.
ORCA는 최적화가 끝나면 같은 구조를 입력이름.xyz로도 저장하니, 출력 파일만 남은 경우에 이 방법을 쓴다.
여러 계산을 CSV 표로
계산마다 디렉터리가 하나씩 있고 그 안에 출력이 있다고 하자. 끝나지 않은 계산은 건너뛰고, 나머지의 에너지를 모은 뒤 가장 낮은 값 기준 상대 에너지(1 Eh = 627.5095 kcal/mol)를 붙인다.
#!/bin/bash
# collect.sh: final energy of every job directory into a CSV
set -euo pipefail
echo "name,energy_Eh" > energies.csv
for out in jobs/*/*.out; do
name=$(basename "$(dirname "$out")")
if ! grep -q 'ORCA TERMINATED NORMALLY' "$out"; then
echo "skip $name: not finished" >&2
continue
fi
e=$(awk '/FINAL SINGLE POINT ENERGY/ {e = $5} END {print e}' "$out")
echo "$name,$e" >> energies.csv
done
awk -F, 'NR > 1 {name[NR] = $1; e[NR] = $2; if (NR == 2 || $2 < min) min = $2}
END {print "name,energy_Eh,rel_kcal"
for (i = 2; i <= NR; i++) printf "%s,%.8f,%.2f\n", name[i], e[i], (e[i] - min) * 627.5095}' \
energies.csv > table.csv$ bash collect.sh
skip 05-killed: not finished
$ column -t -s, table.csv
name energy_Eh rel_kcal
01-single-point -76.32127441 65.95
02-opt-freq -76.42636946 0.00
03-dlpno-ccsdt -76.33195934 59.24
04-solvation-smd -76.33297155 58.61여기 쓴 네 계산은 같은 물 분자를 서로 다른 방법과 기저로 계산한 것이라, 상대 에너지 자체에는 화학적 의미가 없다. 의미가 있는 비교는 같은 B3LYP/def2-SVP에서 기체상(01)과 SMD 용매(04)를 비교한 −7.34 kcal/mol 하나다. 실제로는 같은 방법으로 계산한 이성질체나 배좌 디렉터리에 이 스크립트를 돌린다.
Windows에서 만들거나 옮긴 파일은 줄 끝이 CRLF(\r\n)일 수 있다. 그러면 마지막 열에 \r이 붙어
CSV에 보이지 않는 문자가 섞이고 $ 앵커가 맞지 않는다. cat -A 파일로 줄 끝에 ^M이 보이면
tr -d '\r' < 원본 > 새파일이나 dos2unix로 고친다.
Quantum ESPRESSO
pw.x는 SCF 반복마다 total energy를 쓰고, 수렴한 값에만 줄 맨 앞에 !를 붙인다.
그래서 ^!로 찾으면 중간값이 섞이지 않는다. 아래는
QE 가이드의 실리콘 SCF 예제 출력이다. 단위는 Ry이고 1 Ry = 13.605693 eV다.
$ grep '^!' si.scf.out
! total energy = -93.45273690 Ry
$ grep '^!' si.scf.out | tail -1 | awk '{printf "%.8f Ry = %.6f eV\n", $5, $5 * 13.605693}'
-93.45273690 Ry = -1271.489248 eV
$ grep -c 'convergence has been achieved' si.scf.out
1
$ grep -c 'convergence NOT achieved' si.scf.out
0
$ grep 'JOB DONE' si.scf.out
JOB DONE.
vc-relax는 이온 스텝마다 SCF를 새로 돌리므로 ! 줄이 여러 개다. 압력은 P= 줄의 마지막 열이다.
$ grep '^!' si.vcrelax.out | awk '{printf "%d %.8f\n", NR, $5}'
1 -93.44443849
2 -93.45157833
3 -93.45321788
4 -93.45325068
5 -93.45325084
6 -93.45327391
$ grep 'P=' si.vcrelax.out | awk '{print $NF}'
108.71
43.06
5.66
0.39
0.00
-0.03
$ grep -E 'bfgs converged|Final enthalpy' si.vcrelax.out
bfgs converged in 5 scf cycles and 4 bfgs steps
Final enthalpy = -93.4532508364 Ry
여섯째 값은 BFGS가 끝난 뒤 최종 셀에서 다시 돌린 SCF다. 셀이 바뀌면 평면파 기저도 바뀌므로 값이 조금 다르다.
최적화 결과로는 Final enthalpy를 쓰고, 마지막 SCF와 차이가 크면 컷오프가 부족하다는 신호로 본다.
LAMMPS thermo 표
log.lammps에서 thermo 출력은 Step으로 시작하는 머리글 줄과 Loop time 줄 사이에 있다.
run 명령이 여러 번이면 이 구간도 여러 번 나온다. 아래는
LAMMPS 가이드의 Cu 슬랩 예제 로그로, 평형 10000 스텝과 생산 5000 스텝 두 구간이 있다.
$ awk '$1 == "Step" {on = 1; next} /^Loop time/ {on = 0} on' log.lammps | head -4
0 253.8184 -5749.9163 -5695.3556 -5778.9907
500 253.05769 -5694.4003 -5640.0031 -9546.1221
1000 267.19235 -5695.7562 -5638.3206 1729.171
1500 253.1248 -5692.6357 -5638.2241 -9774.9528
머리글은 앞에 공백이 있어서 /^Step/으로는 걸리지 않는다. 첫 열을 비교하는 $1 == "Step"이 공백에 상관없이 맞는다.
머리글을 한 번만 # 주석으로 남기면 numpy의 loadtxt나 gnuplot이 바로 읽는 파일이 된다.
$ awk '$1 == "Step" {on = 1; if (!hdr++) print "#", $0; next} /^Loop time/ {on = 0} on' log.lammps > thermo.dat
$ head -3 thermo.dat
# Step Temp PotEng TotEng Press
0 253.8184 -5749.9163 -5695.3556 -5778.9907
500 253.05769 -5694.4003 -5640.0031 -9546.1221
두 구간이 이어 붙어서 스텝 번호가 0부터 다시 시작한다(이 입력은 생산 전에 reset_timestep 0을 했다).
생산 구간만 필요하면 머리글을 셀 때 두 번째 구간에서만 켠다.
$ awk '$1 == "Step" {run++; on = (run == 2); next} /^Loop time/ {on = 0} on' log.lammps | head -3
0 247.4071 -5695.3154 -5642.1329 -4114.9433
500 254.87546 -5696.5557 -5641.7678 -4072.5059
1000 242.9701 -5696.8721 -5644.6433 -4897.1177 xtb
xtb는 요약 상자 두 곳에 같은 에너지를 쓴다. 아래는 xtb 6.7.1로 물의 GFN2-xTB 단일점(xtb_sp.out)과
메탄올 구조 최적화(xtb_opt.out)를 계산한 출력이다.
TOTAL ENERGY 줄은 마지막에 한 번만 나오므로 최적화 출력에서도 최종값이다.
$ grep -E ':: total energy|TOTAL ENERGY' xtb_sp.out
:: total energy -5.070222286727 Eh ::
| TOTAL ENERGY -5.070222286727 Eh |
$ awk '/TOTAL ENERGY/ {print $4}' xtb_opt.out
-8.226118421023
$ grep 'GEOMETRY OPTIMIZATION CONVERGED' xtb_opt.out
*** GEOMETRY OPTIMIZATION CONVERGED AFTER 5 ITERATIONS ***
xtb는 normal termination of xtb를 표준 출력이 아니라 표준 오류로 쓴다. xtb mol.xyz > sp.out으로만
저장하면 이 줄이 파일에 없으니, 종료 확인에 쓰려면 > sp.out 2>&1로 함께 저장한다.
Python으로 넘어갈 때
grep과 awk는 값 한두 개를 뽑는 데 알맞다. 여러 블록에 걸친 값(진동 모드 벡터, 원자별 전하, 궤도 에너지)이나
프로그램 버전마다 형식이 바뀌는 부분을 다뤄야 한다면 파서 라이브러리를 쓰는 편이 덜 깨진다.
cclib은 ORCA, Gaussian 등 양자화학 출력을 읽어 numpy 배열로 주고,
ASE의 ase.io.read는 QE 출력과 LAMMPS 덤프를 포함해 구조 파일 대부분을 읽는다.