Line data Source code
1 : /*----------------------------------------------------------------------------*/
2 : /* CP2K: A general program to perform molecular dynamics simulations */
3 : /* Copyright 2000-2026 CP2K developers group <https://cp2k.org> */
4 : /* */
5 : /* SPDX-License-Identifier: GPL-2.0-or-later */
6 : /*----------------------------------------------------------------------------*/
7 :
8 : #include "libcp2k.h"
9 : #include "mpiwrap/cp_mpi.h"
10 : #include <math.h>
11 : #include <stdio.h>
12 : #include <stdlib.h>
13 :
14 : /*******************************************************************************
15 : * \brief Unit test of the C-interface provided via libcp2k.h
16 : * \author Ole Schuett
17 : ******************************************************************************/
18 2 : int main() {
19 :
20 2 : printf("Unit test starts ...\n");
21 :
22 : // test cp2k_get_version()
23 2 : printf("Testing cp_c_get_version(): ");
24 2 : char version_str[100];
25 2 : cp2k_get_version(version_str, 100);
26 2 : printf("%s.\n", version_str);
27 :
28 2 : cp2k_init();
29 2 : const int rank = cp_mpi_comm_rank(cp_mpi_get_comm_world());
30 :
31 : // create simple input file
32 2 : const char *inp_fn = "H2.inp";
33 2 : FILE *f;
34 2 : if (rank == 0) {
35 1 : f = fopen(inp_fn, "w");
36 1 : fprintf(f, "&FORCE_EVAL\n");
37 1 : fprintf(f, " METHOD Quickstep\n");
38 1 : fprintf(f, " STRESS_TENSOR ANALYTICAL\n");
39 1 : fprintf(f, " &DFT\n");
40 1 : fprintf(f, " BASIS_SET_FILE_NAME BASIS_SET\n");
41 1 : fprintf(f, " POTENTIAL_FILE_NAME POTENTIAL\n");
42 1 : fprintf(f, " LSD\n");
43 1 : fprintf(f, " &MGRID\n");
44 1 : fprintf(f, " CUTOFF 140\n");
45 1 : fprintf(f, " &END MGRID\n");
46 1 : fprintf(f, " &QS\n");
47 1 : fprintf(f, " EPS_DEFAULT 1.0E-8\n");
48 1 : fprintf(f, " &END QS\n");
49 1 : fprintf(f, " &SCF\n");
50 1 : fprintf(f, " EPS_DIIS 0.1\n");
51 1 : fprintf(f, " EPS_SCF 1.0E-4\n");
52 1 : fprintf(f, " IGNORE_CONVERGENCE_FAILURE\n");
53 1 : fprintf(f, " MAX_DIIS 4\n");
54 1 : fprintf(f, " MAX_SCF 3\n");
55 1 : fprintf(f, " SCF_GUESS atomic\n");
56 1 : fprintf(f, " &END SCF\n");
57 1 : fprintf(f, " &XC\n");
58 1 : fprintf(f, " &XC_FUNCTIONAL Pade\n");
59 1 : fprintf(f, " &END XC_FUNCTIONAL\n");
60 1 : fprintf(f, " &END XC\n");
61 1 : fprintf(f, " &END DFT\n");
62 1 : fprintf(f, " &SUBSYS\n");
63 1 : fprintf(f, " &CELL\n");
64 1 : fprintf(f, " ABC 8.0 4.0 4.0\n");
65 1 : fprintf(f, " &END CELL\n");
66 1 : fprintf(f, " &COORD\n");
67 1 : fprintf(f, " H 0.000000 0.000000 0.000000\n");
68 1 : fprintf(f, " H 1.000000 0.000000 0.000000\n");
69 1 : fprintf(f, " &END COORD\n");
70 1 : fprintf(f, " &KIND H\n");
71 1 : fprintf(f, " BASIS_SET DZV-GTH-PADE\n");
72 1 : fprintf(f, " POTENTIAL GTH-PADE-q1\n");
73 1 : fprintf(f, " &END KIND\n");
74 1 : fprintf(f, " &END SUBSYS\n");
75 1 : fprintf(f, "&END FORCE_EVAL\n");
76 1 : fprintf(f, "&GLOBAL\n");
77 1 : fprintf(f, " PRINT_LEVEL SILENT\n");
78 1 : fprintf(f, " PROJECT libcp2k_unittest_H2\n");
79 1 : fprintf(f, "&END GLOBAL\n");
80 1 : fclose(f);
81 : }
82 2 : cp_mpi_barrier(cp_mpi_get_comm_world());
83 :
84 : // use input file to create a force environment
85 2 : force_env_t force_env;
86 2 : cp2k_create_force_env(&force_env, inp_fn, "__STD_OUT__");
87 2 : int scf_status = 99;
88 2 : cp2k_get_scf_convergence(force_env, &scf_status);
89 2 : if (scf_status != -1) {
90 0 : printf("SCF status must be unavailable before calculation\n");
91 0 : return (-1);
92 : }
93 2 : cp2k_calc_energy_force(force_env);
94 2 : cp2k_get_scf_convergence(force_env, &scf_status);
95 2 : if (scf_status != 0) {
96 0 : printf("The deliberately truncated SCF must report non-convergence\n");
97 0 : return (-1);
98 : }
99 :
100 : // Stress is a column-major, pressure-positive potential tensor.
101 2 : double stress[9];
102 2 : int stress_available = 0;
103 2 : cp2k_get_stress_tensor(force_env, stress, &stress_available);
104 2 : if (!stress_available) {
105 0 : printf("Missing analytical stress\n");
106 0 : return (-1);
107 : }
108 8 : for (int i = 0; i < 3; ++i) {
109 24 : for (int j = 0; j < 3; ++j) {
110 18 : if (!isfinite(stress[3 * j + i]) ||
111 18 : fabs(stress[3 * j + i] - stress[3 * i + j]) > 1e-10) {
112 0 : printf("Invalid stress tensor\n");
113 0 : return (-1);
114 : }
115 : }
116 : }
117 :
118 : // check energy
119 2 : double energy;
120 2 : cp2k_get_potential_energy(force_env, &energy);
121 2 : printf("\n ENERGY: %.12f\n", energy);
122 2 : if (fabs(-1.118912797546392 - energy) / fabs(energy) > 1e-13) {
123 0 : printf("Wrong energy\n");
124 0 : return (-1);
125 : }
126 :
127 : // Geometry and velocity updates invalidate the last SCF status.
128 2 : double positions[6], cell[9], velocities[6] = {0};
129 2 : cp2k_get_positions(force_env, positions, 6);
130 2 : cp2k_set_positions(force_env, positions, 6);
131 2 : cp2k_get_scf_convergence(force_env, &scf_status);
132 2 : if (scf_status != -1)
133 : return (-1);
134 2 : cp2k_calc_energy(force_env);
135 2 : cp2k_get_cell(force_env, cell);
136 2 : cp2k_set_cell(force_env, cell);
137 2 : cp2k_get_scf_convergence(force_env, &scf_status);
138 2 : if (scf_status != -1)
139 : return (-1);
140 2 : cp2k_calc_energy(force_env);
141 2 : cp2k_set_velocities(force_env, velocities, 6);
142 2 : cp2k_get_scf_convergence(force_env, &scf_status);
143 2 : if (scf_status != -1)
144 : return (-1);
145 2 : cp2k_destroy_force_env(force_env);
146 : // A library caller has not already opened output in the Fortran runtime.
147 : // run_input must create, close, and subsequently append to the named file.
148 2 : const char *run_out = "libcp2k_unittest_run.out";
149 2 : cp2k_run_input_comm(inp_fn, run_out,
150 2 : cp_mpi_comm_c2f(cp_mpi_get_comm_world()));
151 2 : cp_mpi_barrier(cp_mpi_get_comm_world());
152 2 : long first_size = 0;
153 2 : if (rank == 0) {
154 1 : f = fopen(run_out, "r");
155 1 : if (f == NULL) {
156 0 : printf("run_input did not create its output file\n");
157 0 : return (-1);
158 : }
159 1 : fseek(f, 0, SEEK_END);
160 1 : first_size = ftell(f);
161 1 : fclose(f);
162 : }
163 2 : cp2k_run_input(inp_fn, run_out);
164 2 : cp_mpi_barrier(cp_mpi_get_comm_world());
165 2 : if (rank == 0) {
166 1 : f = fopen(run_out, "r");
167 1 : if (f == NULL) {
168 0 : printf("run_input output file disappeared\n");
169 0 : return (-1);
170 : }
171 1 : fseek(f, 0, SEEK_END);
172 1 : const long second_size = ftell(f);
173 1 : fclose(f);
174 1 : if (first_size <= 0 || second_size <= first_size) {
175 0 : printf("run_input did not append output\n");
176 0 : return (-1);
177 : }
178 : }
179 :
180 : // run_input must leave the surrounding library runtime usable.
181 2 : cp2k_create_force_env(&force_env, inp_fn, "__STD_OUT__");
182 2 : cp2k_calc_energy_force(force_env);
183 2 : cp2k_get_potential_energy(force_env, &energy);
184 2 : if (fabs(-1.118912797546392 - energy) / fabs(energy) > 1e-13) {
185 0 : printf("Wrong energy after run_input\n");
186 0 : return (-1);
187 : }
188 2 : cp2k_destroy_force_env(force_env);
189 :
190 : // clean up
191 2 : cp2k_finalize();
192 2 : if (rank == 0) {
193 1 : remove(inp_fn);
194 1 : remove(run_out);
195 : }
196 :
197 2 : printf("Unit test finished, found no errors\n");
198 2 : return (0);
199 : }
200 :
201 : // EOF
|
