Actual source code: ex59.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: static char help[] = "1-D discrete Kohn-Sham model solved as a Nonlinear Eigenvalue Problem with eigenvector dependence (NEPv).\n\n"
12: "Implemented by embedding an EPS inside a nonlinear SNES loop.\n\n"
13: "The Kohn-Sham (KS) equations simplify a complex, interacting many-electron system by mapping it\n"
14: "onto an equivalent system of noninteracting particles. The goal is to find the single-particle\n"
15: "wavefunctions (eigenvectors) and orbital energies (eigenvalues) that self-consistently describe the ground state.\n\n"
16: "The nonlinear Hamiltonian matrix is defined as H(X) = L + alpha * Diag(L^-1 * rho(X)), where:\n"
17: " L = 1-D discrete Laplacian operator matrix.\n"
18: " rho = electronic density vector computed from the current block of eigenvectors X as rho = diag(X * X').\n\n"
19: "The command line options are:\n"
20: " -n <n>, where <n> = number of grid points.\n"
21: " -k <k>, where <k> = number of eigenvectors to compute.\n"
22: " -alpha <alpha>, where <alpha> = real scaling parameter controlling the nonlinearity.\n"
23: " -verbose, to print the converged Hamiltonian diagonal and its sum.\n\n";
25: #include <slepceps.h>
27: /*
28: Context structure to store the necessary data regarding the discrete Kohn-Sham (DKS) problem
29: as well as the objects needed for the nonlinear iteration.
30: */
31: typedef struct {
32: // DKS context
34: // Inputs
35: PetscInt n; // Spatial mesh size
36: PetscInt k; // Number of eigenvectors to compute
37: PetscReal alpha; // Parameter controlling the nonlinearity
39: // Constants
40: Mat L; // 1D discrete Laplacian matrix
42: // Working vectors and base solvers
43: KSP ksp; // Linear solver (to apply L^-1)
44: Vec rho; // Vector to store the electronic density
45: Vec z; // Intermediate vector to store the result z = L^-1 * rho
46: Vec xr; // Temporary vector to extract eigenvectors during each iteration
48: // Nonlinear iteration context
49: EPS eps; // Solver for each H
50: Mat H; // Hamiltonian matrix
51: BV X; // Block of eigenvectors
52: Vec *initial_space; // Array of vectors to store the EPS initial space
53: } DKSContext;
55: /*
56: Initialize the context for the 1D discrete Kohn-Sham problem.
58: This function allocates the necessary memory, builds the 1D Laplacian operator,
59: and configures the linear solver (KSP) to apply L^-1.
61: Arguments:
62: comm - MPI communicator
63: n - Spatial mesh size
64: k - Number of eigenvectors to compute
65: alpha - Parameter controlling the nonlinearity
66: ctx_out - Output pointer where the created context will be stored
67: */
68: PetscErrorCode DKSCreateContext(MPI_Comm comm,PetscInt n,PetscInt k,PetscReal alpha,DKSContext **ctx_out)
69: {
70: DKSContext *ctx;
71: PetscInt Istart,Iend,i;
72: PC pc;
74: PetscFunctionBeginUser;
75: PetscCall(PetscNew(&ctx));
76: ctx->n=n;
77: ctx->k=k;
78: ctx->alpha=alpha;
80: // 1. Create and fill the Laplacian L
81: PetscCall(MatCreate(comm,&ctx->L));
82: PetscCall(MatSetSizes(ctx->L,PETSC_DECIDE,PETSC_DECIDE,n,n));
83: PetscCall(MatSetFromOptions(ctx->L));
85: PetscCall(MatGetOwnershipRange(ctx->L,&Istart,&Iend));
86: for (i=Istart;i<Iend;i++) {
87: PetscCall(MatSetValue(ctx->L,i,i,2.0,INSERT_VALUES));
88: if (i>0) PetscCall(MatSetValue(ctx->L,i,i-1,-1.0,INSERT_VALUES));
89: if (i<n-1) PetscCall(MatSetValue(ctx->L,i,i+1,-1.0,INSERT_VALUES));
90: }
91: PetscCall(MatAssemblyBegin(ctx->L,MAT_FINAL_ASSEMBLY));
92: PetscCall(MatAssemblyEnd(ctx->L,MAT_FINAL_ASSEMBLY));
94: // 2. Configure the KSP
95: PetscCall(KSPCreate(comm,&ctx->ksp));
96: PetscCall(KSPSetOperators(ctx->ksp,ctx->L,ctx->L));
97: PetscCall(KSPSetType(ctx->ksp,KSPPREONLY));
98: PetscCall(KSPGetPC(ctx->ksp,&pc));
99: PetscCall(PCSetType(pc,PCCHOLESKY));
101: if (PetscDefined(HAVE_MUMPS)) PetscCall(PCFactorSetMatSolverType(pc,MATSOLVERMUMPS));
103: PetscCall(KSPSetFromOptions(ctx->ksp));
105: // 3. Create internal work vectors
106: PetscCall(MatCreateVecs(ctx->L,&ctx->rho,&ctx->z));
108: *ctx_out=ctx;
109: PetscFunctionReturn(PETSC_SUCCESS);
110: }
112: /*
113: Configure the SLEPc objects necessary for the nonlinear iteration.
115: This function prepares the Hamiltonian matrix by cloning the structure of the Laplacian,
116: initializes the basis vectors (BV) block for the eigenvectors, and configures the eigenvalue
117: solver (EPS) by defining the problem type and adjusting its tolerance.
119: Arguments:
120: ctx - Pointer to the previously initialized DKS context
121: tol - Desired tolerance for the internal eigenvalue solver
122: */
123: PetscErrorCode DKSSetUpIterationWorkspace(DKSContext *ctx,PetscReal tol)
124: {
125: MPI_Comm comm;
127: PetscFunctionBeginUser;
128: comm = PetscObjectComm((PetscObject)ctx->L);
129: /* Prepare the H matrix */
130: PetscCall(MatDuplicate(ctx->L,MAT_DO_NOT_COPY_VALUES,&ctx->H));
132: /* Configure the eigenvector block X */
133: PetscCall(BVCreate(comm,&ctx->X));
134: PetscCall(BVSetSizesFromVec(ctx->X,ctx->rho,ctx->k));
135: PetscCall(BVSetFromOptions(ctx->X));
137: /* Configure SLEPc EPS */
138: PetscCall(EPSCreate(comm,&ctx->eps));
139: PetscCall(EPSSetProblemType(ctx->eps,EPS_HEP));
140: PetscCall(EPSSetWhichEigenpairs(ctx->eps,EPS_SMALLEST_REAL));
141: PetscCall(EPSSetDimensions(ctx->eps,ctx->k,PETSC_DECIDE,PETSC_DECIDE));
142: PetscCall(EPSSetTolerances(ctx->eps,tol,PETSC_DECIDE));
143: PetscCall(EPSSetFromOptions(ctx->eps));
145: /* Pre-allocate temporary vectors for the nonlinear iteration loop */
146: PetscCall(VecDuplicateVecs(ctx->rho,ctx->k,&ctx->initial_space));
147: PetscCall(MatCreateVecs(ctx->H,&ctx->xr,NULL));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: /*
152: Generate the initial guess X0 using the exact eigenvectors of the 1D Laplacian.
154: This function fills the vector block X with the initial guess, which greatly
155: improves the convergence of the nonlinear iteration loop.
157: Formula used: vv = [1:n]'/(n+1)*pi; X0 = sin(vv * [1:k])*sqrt(2/(n+1));
159: Arguments:
160: ctx - Pointer to the previously initialized DKS context
161: */
162: PetscErrorCode DKSGenerateInitialGuess(DKSContext *ctx)
163: {
164: PetscInt Istart,Iend,i,j;
165: PetscReal h_val,norm_factor,val;
166: Vec col;
167: PetscScalar *x_local;
169: PetscFunctionBeginUser;
170: // Mathematical constants of the formula
171: // X0 = sin( [1:n]' * [1:k] * pi/(n+1) ) * sqrt(2/(n+1));
172: h_val=PETSC_PI/(ctx->n+1.0);
173: norm_factor=PetscSqrtReal(2.0/(ctx->n+1.0));
175: // Iterate over each column (eigenvector)
176: for (j=0;j<ctx->k;j++) {
177: PetscCall(BVGetColumn(ctx->X,j,&col));
178: PetscCall(VecGetOwnershipRange(col,&Istart,&Iend));
179: PetscCall(VecGetArray(col,&x_local));
181: for (i=Istart;i<Iend;i++) {
182: // i and j start at 0 in C, so we add 1 for the mathematical formula
183: val=PetscSinReal((i+1.0)*(j+1.0)*h_val)*norm_factor;
184: x_local[i-Istart]=val;
185: }
187: PetscCall(VecRestoreArray(col,&x_local));
188: PetscCall(BVRestoreColumn(ctx->X,j,&col));
189: }
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: /*
194: Compute the eigenvector-dependent term from the eigenvectors.
196: In the context of the Kohn-Sham model, this function computes the electronic
197: density rho(X) as the sum of the squares of the components of each eigenvector.
198: The formula used is: rho(X) = diag(X * X').
200: Arguments:
201: ctx - Pointer to the initialized DKS context
202: rho_out - Pre-created vector where the computed term will be stored
203: */
204: PetscErrorCode DKSComputeEigvecDependentTerm(DKSContext *ctx,Vec rho_out)
205: {
206: PetscInt i,j,n_loc;
207: Vec col;
208: const PetscScalar *x_local;
209: PetscScalar *rho_local;
211: PetscFunctionBeginUser;
212: PetscCall(VecSet(rho_out,0.0));
213: PetscCall(VecGetLocalSize(rho_out,&n_loc));
214: PetscCall(VecGetArray(rho_out,&rho_local));
216: for (j=0;j<ctx->k;j++) {
217: PetscCall(BVGetColumn(ctx->X,j,&col));
218: PetscCall(VecGetArrayRead(col,&x_local));
219: for (i=0;i<n_loc;i++) {
220: rho_local[i]+=x_local[i]*PetscConj(x_local[i]);
221: }
223: PetscCall(VecRestoreArrayRead(col,&x_local));
224: PetscCall(BVRestoreColumn(ctx->X,j,&col));
225: }
227: PetscCall(VecRestoreArray(rho_out,&rho_local));
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: /*
232: Build the Hamiltonian from an input eigenvector-dependent term.
234: This function solves the linear system L * z = rho_in using the configured
235: KSP solver. Then, it assembles the updated Hamiltonian matrix stored in the context
236: using the formula: H = L + alpha * Diag(z), where z = L^-1 * rho_in.
238: Arguments:
239: ctx - Pointer to the initialized DKS context (ctx->H and ctx->z are updated)
240: rho_in - Vector with the proposed input eigenvector-dependent term
241: */
242: PetscErrorCode DKSBuildHamiltonian(DKSContext *ctx,Vec rho_in)
243: {
244: PetscFunctionBeginUser;
245: // 1. Solve L * z = rho_in
246: PetscCall(KSPSolve(ctx->ksp,rho_in,ctx->z));
248: // 2. Build H = L + alpha * Diag(z)
249: PetscCall(MatCopy(ctx->L,ctx->H,SAME_NONZERO_PATTERN));
250: PetscCall(VecScale(ctx->z,ctx->alpha));
251: PetscCall(MatDiagonalSet(ctx->H,ctx->z,ADD_VALUES));
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: /*
256: Evaluation function (Callback) for the SNES nonlinear solver.
258: In each SNES iteration, this function receives a proposed eigenvector-dependent
259: term (rho_in), builds the Hamiltonian, solves the eigenvalue equation, computes the
260: resulting term, and returns the residual F = rho_out - rho_in.
262: Arguments:
263: snes - The nonlinear solver context
264: rho_in - Input state (eigenvector-dependent term) proposed by SNES
265: F - Vector where the computed residual will be stored
266: ctx_void - Pointer to the user's DKS context (DKSContext)
267: */
268: PetscErrorCode DKSNonlinearIteration(SNES snes,Vec rho_in,Vec F,void *ctx_void)
269: {
270: DKSContext *ctx=(DKSContext*)ctx_void;
271: PetscInt j,nconv;
272: Vec col;
273: MPI_Comm comm;
275: PetscFunctionBeginUser;
276: comm = PetscObjectComm((PetscObject)ctx->L);
277: /* 1. Build the physics (Hamiltonian) with the guess proposed by SNES */
278: PetscCall(DKSBuildHamiltonian(ctx,rho_in));
280: /* 2. Solve with the current H */
281: // Pass our newly built H matrix to SLEPc
282: PetscCall(EPSSetOperators(ctx->eps,ctx->H,NULL));
284: // Extract the eigenvectors from the previous iteration to use as the initial space
285: for (j=0; j<ctx->k; j++) {
286: PetscCall(BVGetColumn(ctx->X,j,&col));
287: PetscCall(VecCopy(col,ctx->initial_space[j]));
288: PetscCall(BVRestoreColumn(ctx->X,j,&col));
289: }
291: // Inject the initial space into the EPS
292: PetscCall(EPSSetInitialSpace(ctx->eps,ctx->k,ctx->initial_space));
294: PetscCall(EPSSolve(ctx->eps));
296: // (Safety check: verify that SLEPc has not failed internally)
297: PetscCall(EPSGetConverged(ctx->eps,&nconv));
298: PetscCheck(nconv>=ctx->k,comm,PETSC_ERR_NOT_CONVERGED,"SLEPc only converged %" PetscInt_FMT " out of %" PetscInt_FMT " eigenvalues in this SNES iteration",nconv,ctx->k);
300: /* 3. Extract the new eigenvectors and store them in ctx->X */
301: for (j=0;j<ctx->k;j++) {
302: PetscCall(EPSGetEigenvector(ctx->eps,j,ctx->xr,NULL)); // Extracts the j-th vector
303: PetscCall(BVGetColumn(ctx->X,j,&col)); // Gets column j of X
304: PetscCall(VecCopy(ctx->xr,col)); // Copies the data
305: PetscCall(BVRestoreColumn(ctx->X,j,&col)); // Restores the column
306: }
308: /* 4. Compute the new eigenvector-dependent term */
309: // We use ctx->rho as a temporary work vector to store rho_out
310: PetscCall(DKSComputeEigvecDependentTerm(ctx,ctx->rho));
312: /* 5. Calculate the nonlinear residual: F = rho_in - rho_out */
313: PetscCall(VecWAXPY(F,-1.0,ctx->rho,rho_in));
314: PetscFunctionReturn(PETSC_SUCCESS);
315: }
317: /*
318: Free the memory associated with the DKS context.
320: This function destroys all internal PETSc objects created during
321: initialization and frees the memory of the main structure.
323: Arguments:
324: ctx - Pointer to the DKS context (set to NULL upon completion)
325: */
326: PetscErrorCode DKSDestroyContext(DKSContext **ctx)
327: {
328: PetscFunctionBeginUser;
329: if (!*ctx) PetscFunctionReturn(PETSC_SUCCESS);
331: PetscCall(MatDestroy(&(*ctx)->L));
332: PetscCall(KSPDestroy(&(*ctx)->ksp));
333: PetscCall(VecDestroy(&(*ctx)->rho));
334: PetscCall(VecDestroy(&(*ctx)->z));
336: PetscCall(MatDestroy(&(*ctx)->H));
337: PetscCall(BVDestroy(&(*ctx)->X));
338: PetscCall(EPSDestroy(&(*ctx)->eps));
340: PetscCall(VecDestroy(&(*ctx)->xr));
341: PetscCall(VecDestroyVecs((*ctx)->k, &(*ctx)->initial_space));
343: PetscCall(PetscFree(*ctx));
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: int main(int argc,char **argv)
348: {
349: DKSContext *ctx;
350: Vec H_diag,rho_guess;
351: PetscScalar diagonal_sum,kr;
352: PetscInt n=5,k=3,maxit=100,its,i;
353: PetscInt nconv;
354: PetscReal alpha=0.5,rtol=SLEPC_DEFAULT_TOL,stol=1e-12;
355: SNES snes;
356: Vec F; // Vector to store the residual (F = rho_out - rho_in)
357: PetscBool verbose=PETSC_FALSE;
358: SNESConvergedReason reason;
359: SNESLineSearch linesearch;
361: PetscFunctionBeginUser;
362: PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
364: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"--- DKS SNES NEPv ---\n"));
366: /* 1. Read options from the command line (if provided by the user) */
367: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
368: PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
369: PetscCall(PetscOptionsGetReal(NULL,NULL,"-alpha",&alpha,NULL));
371: /* 2. Initialization of the DKS context */
372: PetscCall(DKSCreateContext(PETSC_COMM_WORLD,n,k,alpha,&ctx));
374: /* 3. Configure general objects for the nonlinear iteration (Matrices, EPS, eigenvectors) */
375: PetscCall(DKSSetUpIterationWorkspace(ctx,rtol*0.1));
377: /* 4. Generate the initial guess for the eigenvectors and their density */
378: PetscCall(DKSGenerateInitialGuess(ctx));
380: /* 5. Calculate the initial density (rho_guess) from X0 */
381: // Since SNES works with densities, the density is extracted from our X0
382: PetscCall(VecDuplicate(ctx->rho,&rho_guess));
383: PetscCall(DKSComputeEigvecDependentTerm(ctx,rho_guess));
385: /* 6. Prepare the residual vector F by cloning the structure of rho */
386: PetscCall(VecDuplicate(ctx->rho,&F));
388: /* 7. Create and set up the nonlinear solver engine (SNES) */
389: PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes));
390: PetscCall(SNESSetFunction(snes,F,DKSNonlinearIteration,ctx));
391: PetscCall(SNESSetTolerances(snes,PETSC_DETERMINE,rtol,stol,maxit,PETSC_DETERMINE));
393: // Default -> NRICHARDSON with step lambda = 1.0 to simulate basic SCF
394: PetscCall(SNESSetType(snes,SNESNRICHARDSON));
395: PetscCall(SNESGetLineSearch(snes,&linesearch));
396: PetscCall(SNESLineSearchSetType(linesearch,SNESLINESEARCHNONE));
397: PetscCall(SNESLineSearchSetDamping(linesearch,1.0));
399: PetscCall(SNESSetFromOptions(snes));
401: PetscCall(SNESGetTolerances(snes,NULL,&rtol,NULL,&maxit,NULL));
402: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Current parameters: n=%" PetscInt_FMT ", k=%" PetscInt_FMT ", alpha=%g, rtol=%g, maxit=%" PetscInt_FMT "\n",n,k,(double)alpha,(double)rtol,maxit));
404: /* 8. Nonlinear iteration loop */
405: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nStarting nonlinear iterations...\n"));
406: PetscCall(SNESSolve(snes,NULL,rho_guess));
408: PetscCall(SNESGetConvergedReason(snes,&reason));
409: PetscCall(SNESGetIterationNumber(snes,&its));
411: /* 9. Result analysis */
412: if (reason>0) {
413: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"--> CONVERGENCE REACHED in %" PetscInt_FMT " iterations (Reason: %s).\n",its,SNESConvergedReasons[reason]));
415: /* Show the final result */
417: PetscCall(EPSGetConverged(ctx->eps,&nconv));
418: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n--- Eigenvalues (%" PetscInt_FMT " found) ---\n",nconv));
420: for (i=0;i<nconv;i++) {
421: PetscCall(EPSGetEigenvalue(ctx->eps,i,&kr,NULL));
422: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Eigenvalue[%" PetscInt_FMT "] = %10.6f\n",i,(double)PetscRealPart(kr)));
423: }
425: PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
427: if (verbose) {
428: PetscCall(MatCreateVecs(ctx->H,NULL,&H_diag));
429: PetscCall(MatGetDiagonal(ctx->H,H_diag));
431: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Diagonal of the converged Hamiltonian:\n"));
432: PetscCall(VecView(H_diag,PETSC_VIEWER_STDOUT_WORLD));
434: PetscCall(VecSum(H_diag,&diagonal_sum));
435: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Sum of the diagonal: %g\n",(double)PetscRealPart(diagonal_sum)));
436: PetscCall(VecDestroy(&H_diag));
437: }
439: } else {
440: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"--> ERROR: SNES did not converge (Reason: %s)\n",SNESConvergedReasons[reason]));
441: }
443: /* 10. Clean up memory */
444: PetscCall(VecDestroy(&F));
445: PetscCall(VecDestroy(&rho_guess));
446: PetscCall(SNESDestroy(&snes));
447: PetscCall(DKSDestroyContext(&ctx));
449: PetscCall(SlepcFinalize());
450: return 0;
451: }
453: /*TEST
455: testset:
456: filter: sed -e "s/1.364212/1.364211/" -e "s/4.864809/4.864808/" -e "s/1.92185[24]/1.921853/" -e "s/2.931104/2.931103/" -e "s/3.95751[46]/3.957515/" -e "s/rtol=1e-05/rtol=1e-08/" -e "s/rtol=1e-16/rtol=1e-08/" -e "s/[0-9]\{1,\} iterations/8 iterations/" -e "s/CONVERGED_SNORM_RELATIVE/CONVERGED_FNORM_RELATIVE/"
457: output_file: output/ex59_1.out
458: test:
459: suffix: 1
460: test:
461: suffix: 2
462: args: -snes_type anderson -snes_anderson_m 7
463: test:
464: suffix: 3
465: args: -snes_type ngmres -npc_snes_type nrichardson -snes_npc_side right
466: test:
467: suffix: 4
468: args: -snes_type composite -snes_composite_type additiveoptimal -snes_composite_sneses anderson,nrichardson -sub_0_snes_anderson_m 7 -sub_0_snes_anderson_beta 0.4
470: TEST*/