Merge pull request #217 from abhirajsharma/master

ChangeLog

@@ -3,6 +3,14 @@
 -Name
 -changes
 
+
+--------------
+September 05, 2024
+Name: Abhiraj Sharma
+Changes: cyclix/cyclix_tools.c, md.c, orbitalElecDens.c, forces.c
+1. Wraparound bug fixed for cyclix
+2. Density extrapolation fixed for cyclix
+
 --------------
 August 16, 2024
 Name: Abhiraj Sharma

src/cyclix/cyclix_tools.c

@@ -211,6 +211,8 @@ void Cart2nonCart_coord_cyclix(const SPARC_OBJ *pSPARC, double *x, double *y, do
     double x1, x2, x3;
     x1 = sqrt((*x) * (*x) + (*y) * (*y));
     x2 = atan2(*y,*x) - pSPARC->twist * (*z);
+    if (*y < 0.0)
+        x2 += 2*M_PI;
     x3 = (*z);
     *x = x1; *y = x2; *z = x3;
 }

src/forces.c

@@ -112,24 +112,6 @@ void Calculate_EGS_Forces(SPARC_OBJ *pSPARC)
 
     // Apply constraint on motion of atoms for Relaxation\MD
     if(pSPARC->RelaxFlag == 1 || pSPARC->MDFlag == 1) {
-        if (pSPARC->CyclixFlag) {
-            double fx, fy;
-            double ty, tz;
-            if (pSPARC->elecgs_Count > 0) {
-                for(i = 0; i < pSPARC->n_atom; i++){
-                    Cart2nonCart_coord(pSPARC, &pSPARC->atom_pos_nm[3*i], &pSPARC->atom_pos_nm[3*i+1], &pSPARC->atom_pos_nm[3*i+2]);
-
-                    ty = (pSPARC->atom_pos_nm[3*i+1] - pSPARC->atom_pos[3*i+1])/pSPARC->range_y;
-                    tz = (pSPARC->atom_pos_nm[3*i+2] - pSPARC->atom_pos[3*i+2])/pSPARC->range_z;
-                    RotMat_cyclix(pSPARC, ty, tz);
-                    fx = pSPARC->forces[3*i]; fy = pSPARC->forces[3*i+1];
-                    pSPARC->forces[3*i] = pSPARC->RotM_cyclix[0] * fx + pSPARC->RotM_cyclix[1] * fy;
-                    pSPARC->forces[3*i+1] = pSPARC->RotM_cyclix[3] * fx + pSPARC->RotM_cyclix[4] * fy;
-
-                    nonCart2Cart_coord(pSPARC, &pSPARC->atom_pos_nm[3*i], &pSPARC->atom_pos_nm[3*i+1], &pSPARC->atom_pos_nm[3*i+2]);
-                }
-            }
-        }
         for(i = 0; i < 3*pSPARC->n_atom; i++)
             pSPARC->forces[i] *= pSPARC->mvAtmConstraint[i];
     }
@@ -1419,7 +1401,7 @@ void Symmetrize_forces(SPARC_OBJ *pSPARC)
     shift_fy /= n_atom;
     shift_fz /= n_atom;
 
-    if (pSPARC->CyclixFlag) {
+    if (pSPARC->CyclixFlag && pSPARC->range_y < 2.0*M_PI) {
         // Do not symmetrize forces in x-y plane for cyclix systems (TODO: But do symmetrize if all the atoms are taken into account)
         shift_fx = 0.0;
         shift_fy = 0.0;

src/include/md.h

@@ -73,10 +73,20 @@ void Calc_vel1_G(SPARC_OBJ *pSPARC);
 void Calc_vel2_G(SPARC_OBJ *pSPARC); 
 
 /*
-* @ brief: function to wrap around atom positions that lie outside main domain in case of PBC and check if the atoms are too close to the boundary in case of bounded domain
+* @ brief: function to check if the atoms are too close to the boundary in case of bounded domain or to each other in general
 */
 void Check_atomlocation(SPARC_OBJ *pSPARC);
 
+/*
+* @ brief: function to wraparound atom positions for PBC
+*/
+void wraparound_dynamics(SPARC_OBJ *pSPARC, double *coord, int opt);
+
+/*
+* @ brief: function to wraparound velocities in MD and displacement vectors in relaxation for PBC
+*/
+void wraparound_velocity(SPARC_OBJ *pSPARC, double shift, int dir, int loc);
+
 /*
  @ brief: function to write all relevant DFT quantities generated during MD simulation
 */

src/initialization.c

@@ -3470,7 +3470,7 @@ void write_output_init(SPARC_OBJ *pSPARC) {
     }
 
     fprintf(output_fp,"***************************************************************************\n");
-    fprintf(output_fp,"*                       SPARC (version August 16, 2024)                      *\n");
+    fprintf(output_fp,"*                       SPARC (version September 05, 2024)                      *\n");
     fprintf(output_fp,"*   Copyright (c) 2020 Material Physics & Mechanics Group, Georgia Tech   *\n");
     fprintf(output_fp,"*           Distributed under GNU General Public License 3 (GPL)          *\n");
     fprintf(output_fp,"*                   Start time: %s                  *\n",c_time_str);

src/makefile

@@ -33,7 +33,7 @@ ENABLE_SIMD_COMPLEX = 1
 # BLASROOT = /usr/lib64
 # LDFLAGS += -L$(BLASROOT)
 
-CPPFLAGS = -Iinclude/ -Ixc/exx/include -Ixc/vdW/d3/include -Ixc/vdW/vdWDF/include -Ixc/mgga/include -IhighT/include -Icyclix/include -Icyclix/cylix_mlff -IhighT/mlff_highT -Imlff/
+CPPFLAGS = -Iinclude/ -Ixc/exx/include -Ixc/vdW/d3/include -Ixc/vdW/vdWDF/include -Ixc/mgga/include -IhighT/include -Icyclix/include -Icyclix/cyclix_mlff -IhighT/mlff_highT -Imlff/
 LDLIBS   = -lrt
 
 ifeq ($(USE_MKL), 1)
@@ -102,7 +102,7 @@ OBJSC = main.o initialization.o readfiles.o atomdata.o parallelization.o relax.o
         cyclix/cyclix_stress.o \
         mlff/covariance_matrix.o mlff/regression.o mlff/spherical_harmonics.o mlff/soap_descriptor.o mlff/sparc_mlff_interface.o \
         mlff/bessel_NR.o mlff/nrutils.o mlff/sparsification.o mlff/tools_mlff.o mlff/mlff_read_write.o \
-        mlff/descriptor.o mlff/ddbp_tools.o highT/mlff_highT/internal_energy_model.o cyclix/cylix_mlff/cyclix_mlff_tools.o \
+        mlff/descriptor.o mlff/ddbp_tools.o highT/mlff_highT/internal_energy_model.o cyclix/cyclix_mlff/cyclix_mlff_tools.o \
         mlff/hnl_soap.o
 
 LIBBASE = ../lib/sparc

src/md.c

@@ -1768,7 +1768,7 @@ void Cart2nonCart(double *gradT, double *carCoord, double *nonCarCoord) {
 }
 
 /*
- @ brief: function to wrap around atom positions that lie outside main domain in case of PBC and check if the atoms are too close to the boundary in case of bounded domain
+* @ brief: function to check if the atoms are too close to the boundary in case of bounded domain or to each other in general
 */
 void Check_atomlocation(SPARC_OBJ *pSPARC) {
     int rank, ityp, i, atm, atm2, count, dir = 0, maxdir = 3, BC;
@@ -1813,10 +1813,22 @@ void Check_atomlocation(SPARC_OBJ *pSPARC) {
 	}
 	free(rc);
 
+	wraparound_dynamics(pSPARC, pSPARC->atom_pos, 1);
+}
+
+/*
+* @ brief: function to wraparound atom positions for PBC
+*/
+void wraparound_dynamics(SPARC_OBJ *pSPARC, double *coord, int opt) {
+
+	int rank, atm, dir = 0, maxdir = 3, BC;
+	double length, coord_temp;// Change tol according to the situation
+	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
+
 	// Convert Cart to nonCart coordinates for non orthogonal cell
     if(pSPARC->cell_typ != 0){
         for(atm = 0; atm < pSPARC->n_atom; atm++)
-	        Cart2nonCart_coord(pSPARC, &pSPARC->atom_pos[3*atm], &pSPARC->atom_pos[3*atm+1], &pSPARC->atom_pos[3*atm+2]);
+	        Cart2nonCart_coord(pSPARC, coord+3*atm, coord+3*atm+1, coord+3*atm+2);
     }
 
 	while(dir < maxdir){
@@ -1834,27 +1846,64 @@ void Check_atomlocation(SPARC_OBJ *pSPARC) {
 		}
 		if(BC == 1){
 			if (!((pSPARC->CyclixFlag) && (dir == 0))) { // if it is in cyclix coordinate and in x (radial) direction, we will not check it
-				for(atm = 0; atm < pSPARC->n_atom; atm++){
-					if(pSPARC->atom_pos[atm * 3 + dir] >= length || pSPARC->atom_pos[atm * 3 + dir] < 0){
-						if(!rank)
-							printf("ERROR: Atom number %d has crossed the boundary in %d direction",atm, dir);
-						exit(EXIT_FAILURE);
+                for(atm = 0; atm < pSPARC->n_atom; atm++){
+                    if(coord[atm * 3 + dir] >= length || coord[atm * 3 + dir] < 0){
+                        if(!rank)
+                            printf("ERROR: Atom number %d has crossed the boundary in %d direction",atm, dir);
+                        exit(EXIT_FAILURE);
+                    }
+                }
+            }
+		} else if(BC == 0){
+			for(atm = 0; atm < pSPARC->n_atom; atm++){
+				coord_temp = *(coord+3*atm+dir);
+            	if (coord_temp < 0.0 || coord_temp >= length)
+            	{
+                	coord_temp = fmod(coord_temp,length);
+					double new_coord = coord_temp + (coord_temp<0.0)*length;
+					double shift = new_coord - *(coord+3*atm+dir);
+					*(coord+3*atm+dir) = new_coord;
+					if (pSPARC->CyclixFlag && opt == 1){
+						wraparound_velocity(pSPARC, shift, dir, 3*atm);
 					}
 				}
 			}
-		}else if(BC == 0){
-		// TODO: Assumed that atom moves less than one cell length, generalize in future
-			for(atm = 0; atm < pSPARC->n_atom; atm++){
-				if(pSPARC->atom_pos[atm * 3 + dir] >= length)
-					pSPARC->atom_pos[atm * 3 + dir] -= length;
-				else if(pSPARC->atom_pos[atm * 3 + dir] < 0)
-					pSPARC->atom_pos[atm * 3 + dir] += length;
-			}
 		}
 		dir ++;
 	}
 }
 
+/*
+* @ brief: function to wraparound velocities in MD and displacement vectors in relaxation for PBC
+*/
+void wraparound_velocity(SPARC_OBJ *pSPARC, double shift, int dir, int loc) {
+
+	if (dir == 1){
+		if (pSPARC->MDFlag == 1){
+			double vx = pSPARC->ion_vel[loc]; double vy = pSPARC->ion_vel[loc+1];
+			pSPARC->ion_vel[loc] = cos(shift) * vx -sin(shift) * vy;
+			pSPARC->ion_vel[loc+1] = sin(shift) * vx + cos(shift) * vy;
+		}
+		else if (pSPARC->RelaxFlag == 1){
+			double vx = pSPARC->d[loc]; double vy = pSPARC->d[loc+1];
+			pSPARC->d[loc] = cos(shift) * vx -sin(shift) * vy;
+			pSPARC->d[loc+1] = sin(shift) * vx + cos(shift) * vy;
+		}	
+	} else if (dir == 2){
+		if (pSPARC->MDFlag == 1){
+			double vx = pSPARC->ion_vel[loc]; double vy = pSPARC->ion_vel[loc+1];
+			pSPARC->ion_vel[loc] = cos(pSPARC->twist*shift) * vx -sin(pSPARC->twist*shift) * vy;
+			pSPARC->ion_vel[loc+1] = sin(pSPARC->twist*shift) * vx + cos(pSPARC->twist*shift) * vy;
+		}
+		else if (pSPARC->RelaxFlag == 1){
+			double vx = pSPARC->d[loc]; double vy = pSPARC->d[loc+1];
+			pSPARC->d[loc] = cos(pSPARC->twist*shift) * vx -sin(pSPARC->twist*shift) * vy;
+			pSPARC->d[loc+1] = sin(pSPARC->twist*shift) * vx + cos(pSPARC->twist*shift) * vy;
+		}	
+	}
+
+}
+
 /*
  @ brief: function to write all relevant DFT quantities generated during MD simulation
 */

src/orbitalElecDensInit.c

@@ -28,6 +28,8 @@
 #include "electronDensity.h"
 #include "parallelization.h"
 #include "readfiles.h"
+#include "initialization.h"
+#include "md.h"
 #define max(x,y) ((x)>(y)?(x):(y))
 
 
@@ -205,47 +207,93 @@ void Init_electronDensity(SPARC_OBJ *pSPARC) {
  */
 // TODO: Check if the FtF matrix is singular and if it is then decide on how to do extrapolation or not to do at all
 void elecDensExtrapolation(SPARC_OBJ *pSPARC) {
-	// processors that are not in the dmcomm_phi will remain idle
+    // processors that are not in the dmcomm_phi will remain idle
     if (pSPARC->dmcomm_phi == MPI_COMM_NULL) {
         return; 
     }
     int nd, ii, matrank, count = 0, atm;
-    double alpha, beta;
+    double alpha, beta, *coord_temp1, *coord_temp2;
     for (nd = 0; nd < pSPARC->Nd_d; nd++){
         pSPARC->delectronDens_2dt[nd] = pSPARC->delectronDens_1dt[nd];
         pSPARC->delectronDens_1dt[nd] = pSPARC->delectronDens_0dt[nd];
-        pSPARC->delectronDens_0dt[nd] = pSPARC->electronDens[nd] - pSPARC->electronDens_at[nd];
+        double *electronDens = pSPARC->electronDens;        
+        pSPARC->delectronDens_0dt[nd] = electronDens[nd] - pSPARC->electronDens_at[nd];
     }
     if(pSPARC->MDFlag == 1){
-	    for(atm = 0; atm < pSPARC->n_atom; atm++){
-	    	if(pSPARC->MDCount == 1){
-	    		pSPARC->atom_pos_nm[count * 3] = pSPARC->atom_pos[count * 3];
-		    	pSPARC->atom_pos_nm[count * 3 + 1] = pSPARC->atom_pos[count * 3 + 1];
-			    pSPARC->atom_pos_nm[count * 3 + 2] = pSPARC->atom_pos[count * 3 + 2];
-		    	count ++;
-		    } else{
-			    pSPARC->atom_pos_nm[count * 3] += pSPARC->MD_dt * pSPARC->ion_vel[count * 3];
-			    pSPARC->atom_pos_nm[count * 3 + 1] += pSPARC->MD_dt * pSPARC->ion_vel[count * 3 + 1];
-			    pSPARC->atom_pos_nm[count * 3 + 2] += pSPARC->MD_dt * pSPARC->ion_vel[count * 3 + 2];
-			    count ++;
-		    }         
-	    }
+        if(pSPARC->MDCount == 1) {
+            for(atm = 0; atm < pSPARC->n_atom; atm++){
+                pSPARC->atom_pos_nm[count * 3] = pSPARC->atom_pos[count * 3];
+                pSPARC->atom_pos_nm[count * 3 + 1] = pSPARC->atom_pos[count * 3 + 1];
+                pSPARC->atom_pos_nm[count * 3 + 2] = pSPARC->atom_pos[count * 3 + 2];
+                count ++;
+            }   
+        } else{
+            coord_temp1 = (double *) malloc(3*pSPARC->n_atom*sizeof(double));
+            coord_temp2 = (double *) malloc(3*pSPARC->n_atom*sizeof(double));
+
+            for(atm = 0; atm < 3*pSPARC->n_atom; atm++){
+                coord_temp1[atm] = pSPARC->atom_pos_nm[atm];
+                coord_temp2[atm] = pSPARC->atom_pos[atm];
+            }
+            wraparound_dynamics(pSPARC, coord_temp1, 0);
+            if(pSPARC->cell_typ != 0){
+                for(atm = 0; atm < pSPARC->n_atom; atm++){
+                    Cart2nonCart_coord(pSPARC, coord_temp2+3*atm, coord_temp2+3*atm+1, coord_temp2+3*atm+2);
+                    Cart2nonCart_coord(pSPARC, pSPARC->atom_pos_nm+3*atm, pSPARC->atom_pos_nm+3*atm+1, pSPARC->atom_pos_nm+3*atm+2);
+                }
+            }
+            for(atm = 0; atm < 3*pSPARC->n_atom; atm++){
+                pSPARC->atom_pos_nm[atm] += coord_temp2[atm] - coord_temp1[atm];
+            }
+
+
+
+            if(pSPARC->cell_typ != 0){
+                for(atm = 0; atm < pSPARC->n_atom; atm++){
+                    nonCart2Cart_coord(pSPARC, pSPARC->atom_pos_nm+3*atm, pSPARC->atom_pos_nm+3*atm+1, pSPARC->atom_pos_nm+3*atm+2);
+                }
+            }
+
+            free(coord_temp1);
+            free(coord_temp2);
+        }   
     } else if(pSPARC->RelaxFlag == 1){
-        for(atm = 0; atm < pSPARC->n_atom; atm++){
-		    if((pSPARC->elecgs_Count - pSPARC->StressCount) == 1){
-			    pSPARC->atom_pos_nm[count * 3] = pSPARC->atom_pos[count * 3];
-			    pSPARC->atom_pos_nm[count * 3 + 1] = pSPARC->atom_pos[count * 3 + 1];
-			    pSPARC->atom_pos_nm[count * 3 + 2] = pSPARC->atom_pos[count * 3 + 2];
+        if((pSPARC->elecgs_Count - pSPARC->StressCount) == 1) {
+            for(atm = 0; atm < pSPARC->n_atom; atm++){
+                pSPARC->atom_pos_nm[count * 3] = pSPARC->atom_pos[count * 3];
+                pSPARC->atom_pos_nm[count * 3 + 1] = pSPARC->atom_pos[count * 3 + 1];
+                pSPARC->atom_pos_nm[count * 3 + 2] = pSPARC->atom_pos[count * 3 + 2];
                 count ++;
-		    } else{
-			    pSPARC->atom_pos_nm[count * 3] += pSPARC->Relax_fac * pSPARC->d[count * 3] * pSPARC->mvAtmConstraint[count * 3];
-			    pSPARC->atom_pos_nm[count * 3 + 1] += pSPARC->Relax_fac * pSPARC->d[count * 3 + 1] * pSPARC->mvAtmConstraint[count * 3 + 1];
-			    pSPARC->atom_pos_nm[count * 3 + 2] += pSPARC->Relax_fac * pSPARC->d[count * 3 + 2] * pSPARC->mvAtmConstraint[count * 3 + 2];
-			    count ++;
-		    }
-	    }
+            }
+        } else{
+            coord_temp1 = (double *) malloc(3*pSPARC->n_atom*sizeof(double));
+            coord_temp2 = (double *) malloc(3*pSPARC->n_atom*sizeof(double));
+
+            for(atm = 0; atm < 3*pSPARC->n_atom; atm++){
+                coord_temp1[atm] = pSPARC->atom_pos_nm[atm];
+                coord_temp2[atm] = pSPARC->atom_pos[atm];
+            }
+            wraparound_dynamics(pSPARC, coord_temp1, 0);
+            if(pSPARC->cell_typ != 0){
+                for(atm = 0; atm < pSPARC->n_atom; atm++){
+                    Cart2nonCart_coord(pSPARC, coord_temp2+3*atm, coord_temp2+3*atm+1, coord_temp2+3*atm+2);
+                    Cart2nonCart_coord(pSPARC, pSPARC->atom_pos_nm+3*atm, pSPARC->atom_pos_nm+3*atm+1, pSPARC->atom_pos_nm+3*atm+2);
+                }
+            }
+            for(atm = 0; atm < 3*pSPARC->n_atom; atm++){
+                pSPARC->atom_pos_nm[atm] += coord_temp2[atm] - coord_temp1[atm];
+            }
+
+            if(pSPARC->cell_typ != 0){
+                for(atm = 0; atm < pSPARC->n_atom; atm++)
+                    nonCart2Cart_coord(pSPARC, pSPARC->atom_pos_nm+3*atm, pSPARC->atom_pos_nm+3*atm+1, pSPARC->atom_pos_nm+3*atm+2);
+            }
+
+            free(coord_temp1);
+            free(coord_temp2);
+        }   
     }
-	if((pSPARC->elecgs_Count - pSPARC->StressCount) >= 3){ 
+    if((pSPARC->elecgs_Count - pSPARC->StressCount) >= 3){ 
        double *FtF, *Ftf, *s, temp1, temp2, temp3; // 2x2 matrix and 2x1 vectors in (FtF)*(svec) = Ftf
         FtF = (double *)calloc( 4 , sizeof(double) );
         Ftf = (double *)calloc( 2 , sizeof(double) );
@@ -279,7 +327,6 @@ void elecDensExtrapolation(SPARC_OBJ *pSPARC) {
     }
 }
 
-
 /**
  * @brief   initialize Kohn-Sham orbitals.
  */