diff --git a/PuReMD/src/forces.c b/PuReMD/src/forces.c
index bd9fa99da8d4ceb1db3efb0eba17f6140f1690c1..d5132fc003df3326f65ef35b38ef7a8841e2e7ba 100644
--- a/PuReMD/src/forces.c
+++ b/PuReMD/src/forces.c
@@ -347,6 +347,7 @@ void Init_Forces( reax_system *system, control_params *control,
int start_j, end_j;
int btop_j;
#if defined(NEUTRAL_TERRITORY)
+ int mark[6];
int total_cnt[6];
int bin[6];
int total_sum[6];
@@ -373,11 +374,9 @@ void Init_Forces( reax_system *system, control_params *control,
num_hbonds = 0;
btop_i = 0;
renbr = (data->step - data->prev_steps) % control->reneighbor == 0;
-#if defined(NEUTRAL_TERRITORY)
- nt_flag = 1;
-#endif
#if defined(NEUTRAL_TERRITORY)
+ nt_flag = 1;
if( renbr )
{
for ( i = 0; i < 6; ++i )
@@ -410,20 +409,9 @@ void Init_Forces( reax_system *system, control_params *control,
}
H->NT = total_sum[5] + total_cnt[5];
}
-#endif
- /*int local_to_local = 0;
- FILE *fp;*/
-#if defined(NEUTRAL_TERRITORY)
- /*int local_to_nt = 0;
- int nt_to_local = 0;
- int nt_to_nt = 0;
- fp = fopen( "NT_Pairs.txt", "w" );*/
- int mark[6];
mark[0] = mark[1] = 1;
mark[2] = mark[3] = mark[4] = mark[5] = 2;
-#else
- //fp = fopen( "Jacobi_Pairs.txt", "w" );
#endif
for ( i = 0; i < system->N; ++i )
@@ -465,15 +453,12 @@ void Init_Forces( reax_system *system, control_params *control,
ihb = -1;
ihb_top = -1;
- //TODO: undo
if ( local == 1 )
{
H->start[i] = Htop;
H->entries[Htop].j = i;
H->entries[Htop].val = sbp_i->eta;
++Htop;
- /*fprintf( fp, "%d %d\n", atom_i->orig_id, atom_i->orig_id );
- ++local_to_local;*/
if ( control->hbond_cut > 0 )
{
@@ -537,39 +522,35 @@ void Init_Forces( reax_system *system, control_params *control,
if ( local == 1 )
{
/* H matrix entry */
- if ( (far_nbrs->format == HALF_LIST
- && (j < system->n || atom_i->orig_id < atom_j->orig_id))
- || far_nbrs->format == FULL_LIST )
- {
#if defined(NEUTRAL_TERRITORY)
- if( atom_j->nt_dir > 0 || j < system->n )
+ if( atom_j->nt_dir > 0 || (j < system->n && (H->format == SYM_FULL_MATRIX
+ || (H->format == SYM_HALF_MATRIX && i < j))) )
+ {
+ if( j < system->n )
+ {
+ H->entries[Htop].j = j;
+ }
+ else
+ {
+ H->entries[Htop].j = atom_j->pos;
+ }
+
+ if ( control->tabulate == 0 )
+ {
+ H->entries[Htop].val = Compute_H(r_ij, twbp->gamma, workspace->Tap);
+ }
+ else
{
- //fprintf( fp, "%d %d\n", atom_i->orig_id, atom_j->orig_id );
- if( j < system->n )
- {
- //++local_to_local;
- H->entries[Htop].j = j;
- }
- else
- {
- //++local_to_nt;
- H->entries[Htop].j = atom_j->pos;
- }
-
- if ( control->tabulate == 0 )
- {
- H->entries[Htop].val = Compute_H(r_ij, twbp->gamma, workspace->Tap);
- }
- else
- {
- H->entries[Htop].val = Compute_tabH(r_ij, type_i, type_j);
- }
- ++Htop;
+ H->entries[Htop].val = Compute_tabH(r_ij, type_i, type_j);
}
+ ++Htop;
+ }
#else
- /*fprintf( fp, "%d %d\n", atom_i->orig_id, atom_j->orig_id );
+ if ( (far_nbrs->format == HALF_LIST
+ && (j < system->n || atom_i->orig_id < atom_j->orig_id))
+ || far_nbrs->format == FULL_LIST )
+ {
if( j < system->n )
- ++local_to_local;*/
H->entries[Htop].j = j;
if ( control->tabulate == 0 )
{
@@ -580,8 +561,8 @@ void Init_Forces( reax_system *system, control_params *control,
H->entries[Htop].val = Compute_tabH(r_ij, type_i, type_j);
}
++Htop;
-#endif
}
+#endif
/* hydrogen bond lists */
if ( control->hbond_cut > 0 && (ihb == 1 || ihb == 2) &&
@@ -614,9 +595,10 @@ void Init_Forces( reax_system *system, control_params *control,
else if ( local == 2 )
{
/* H matrix entry */
- if( ( atom_j->nt_dir != -1 && mark[atom_i->nt_dir] != mark[atom_j->nt_dir] ) || ( j < system->n && atom_i->nt_dir != 0 ))
+ if( ( atom_j->nt_dir != -1 && mark[atom_i->nt_dir] != mark[atom_j->nt_dir]
+ && ( H->format == SYM_FULL_MATRIX || (H->format == SYM_HALF_MATRIX && atom_i->pos < atom_j->pos)))
+ || ( j < system->n && atom_i->nt_dir != 0 && H->format == SYM_FULL_MATRIX ))
{
- //fprintf( fp, "%d %d\n", atom_i->orig_id, atom_j->orig_id );
if( !nt_flag )
{
nt_flag = 1;
@@ -624,12 +606,10 @@ void Init_Forces( reax_system *system, control_params *control,
}
if( j < system->n )
{
- //++nt_to_local;
H->entries[Htop].j = j;
}
else
{
- //++nt_to_nt;
H->entries[Htop].j = atom_j->pos;
}
@@ -720,30 +700,6 @@ void Init_Forces( reax_system *system, control_params *control,
}
}
-/* int Htot, localtot;
-
- MPI_Allreduce(&Htop, &Htot, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
- MPI_Allreduce(&local_to_local, &localtot, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-
-#if defined(NEUTRAL_TERRITORY)
- int tot2, tot3, tot4;
- MPI_Allreduce(&local_to_nt, &tot2, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
- MPI_Allreduce(&nt_to_local, &tot3, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
- MPI_Allreduce(&nt_to_nt, &tot4, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-#endif
-
- if(system->my_rank == MASTER_NODE )
- {
- fprintf( stdout, "total number of entries across all matrices: %d\n", Htot );
- fprintf( stdout, "total number of local to local entries across all matrices: %d\n", localtot );
-#if defined(NEUTRAL_TERRITORY)
- fprintf( stdout, "total number of local to nt entries across all matrices: %d\n", tot2 );
- fprintf( stdout, "total number of nt to local entries across all matrices: %d\n", tot3 );
- fprintf( stdout, "total number of nt to nt entries across all matrices: %d\n", tot4 );
-#endif
- fflush( stdout );
- }*/
-
workspace->realloc.Htop = Htop;
workspace->realloc.num_bonds = num_bonds;
workspace->realloc.num_hbonds = num_hbonds;
@@ -770,7 +726,6 @@ void Init_Forces( reax_system *system, control_params *control,
Validate_Lists( system, workspace, lists, data->step,
system->n, system->N, system->numH, comm );
- //fclose( fp );
}
@@ -1008,7 +963,8 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
void Estimate_Storages( reax_system *system, control_params *control,
reax_list **lists, int *Htop, int *hb_top,
- int *bond_top, int *num_3body, MPI_Comm comm, int *matrix_dim )
+ int *bond_top, int *num_3body, MPI_Comm comm,
+ int *matrix_dim, int cm_format )
{
int i, j, pj;
int start_i, end_i;
@@ -1032,6 +988,12 @@ void Estimate_Storages( reax_system *system, control_params *control,
memset( bond_top, 0, sizeof(int) * system->total_cap );
*num_3body = 0;
+#if defined(NEUTRAL_TERRITORY)
+ int mark[6];
+ mark[0] = mark[1] = 1;
+ mark[2] = mark[3] = mark[4] = mark[5] = 2;
+#endif
+
for ( i = 0; i < system->N; ++i )
{
atom_i = &(system->my_atoms[i]);
@@ -1053,8 +1015,6 @@ void Estimate_Storages( reax_system *system, control_params *control,
{
local = 2;
cutoff = control->nonb_cut;
- // TODO: Diag entries for NT atoms?
- //++(*Htop);
++(*matrix_dim);
ihb = -1;
}
@@ -1070,7 +1030,7 @@ void Estimate_Storages( reax_system *system, control_params *control,
{
nbr_pj = &( far_nbrs->far_nbr_list[pj] );
j = nbr_pj->nbr;
- if ( far_nbrs->format == HALF_LIST )
+ //if ( far_nbrs->format == HALF_LIST )
{
atom_j = &(system->my_atoms[j]);
}
@@ -1084,17 +1044,19 @@ void Estimate_Storages( reax_system *system, control_params *control,
if ( local == 1 )
{
+#if defined(NEUTRAL_TERRITORY)
+ if( atom_j->nt_dir > 0 || j < system->n )
+ {
+ ++(*Htop);
+ }
+#else
if ( (far_nbrs->format == HALF_LIST
- && (j < system->n || atom_i->orig_id < atom_j->orig_id))
- || far_nbrs->format == FULL_LIST )
+ && (j < system->n || atom_i->orig_id < atom_j->orig_id))
+ || far_nbrs->format == FULL_LIST )
{
-#if defined(NEUTRAL_TERRITORY)
- if( j < system->n || (system->my_atoms[j]).nt_dir != -1 )
-#endif
- {
- ++(*Htop);
- }
+ ++(*Htop);
}
+#endif
/* hydrogen bond lists */
if ( control->hbond_cut > 0.1 && (ihb == 1 || ihb == 2) &&
@@ -1117,15 +1079,13 @@ void Estimate_Storages( reax_system *system, control_params *control,
#if defined(NEUTRAL_TERRITORY)
else if ( local == 2 )
{
-
-#if defined(HALF_LIST)
- if ( j < system->n || atom_i->orig_id < atom_j->orig_id ) //tryQEq ||1
-#endif
+ /*if( ( atom_j->nt_dir != -1 && mark[atom_i->nt_dir] != mark[atom_j->nt_dir]
+ && ( cm_format == SYM_FULL_MATRIX || (cm_format == SYM_HALF_MATRIX && atom_i->pos < atom_j->pos)))
+ || ( j < system->n && atom_i->nt_dir != 0 && cm_format == SYM_FULL_MATRIX ))*/
+ if( ( atom_j->nt_dir != -1 && mark[atom_i->nt_dir] != mark[atom_j->nt_dir] )
+ || ( j < system->n && atom_i->nt_dir != 0 ))
{
- if( j < system->n || (system->my_atoms[j]).nt_dir != -1 )
- {
- ++(*Htop);
- }
+ ++(*Htop);
}
}
#endif
@@ -1172,6 +1132,15 @@ void Estimate_Storages( reax_system *system, control_params *control,
}
}
+#if defined(NEUTRAL_TERRITORY)
+
+ /* Since we don't know the NT atoms' position yet, Htop cannot be calculated accurately
+ Therefore, we assume it is full and divide 2 if necessary */
+ if( cm_format == SYM_HALF_MATRIX )
+ {
+ *Htop = (*Htop + system->n)/2;
+ }
+#endif
*Htop = (int)(MAX( *Htop * SAFE_ZONE, MIN_CAP * MIN_HENTRIES ));
*matrix_dim = (int)(MAX( *matrix_dim * SAFE_ZONE, MIN_CAP ));
for ( i = 0; i < system->n; ++i )
diff --git a/PuReMD/src/forces.h b/PuReMD/src/forces.h
index 10b34395b3af6940a1779255fb8a3ebea7b45976..105f35941c8aa0d791f12b55a59649ba2d4ead8b 100644
--- a/PuReMD/src/forces.h
+++ b/PuReMD/src/forces.h
@@ -31,5 +31,5 @@ void Init_Force_Functions( control_params* );
void Compute_Forces( reax_system*, control_params*, simulation_data*,
storage*, reax_list**, output_controls*, mpi_datatypes* );
void Estimate_Storages( reax_system*, control_params*, reax_list**,
- int*, int*, int*, int*, MPI_Comm, int* );
+ int*, int*, int*, int*, MPI_Comm, int*, int );
#endif
diff --git a/PuReMD/src/init_md.c b/PuReMD/src/init_md.c
index 4e3000b58763a22ae047c22d5806a6f1c903be1a..b8bf0a7238ed0eb91e293cbecebc3d77b1085cc4 100644
--- a/PuReMD/src/init_md.c
+++ b/PuReMD/src/init_md.c
@@ -598,7 +598,7 @@ int Init_Lists( reax_system *system, control_params *control,
{
#if defined(NEUTRAL_TERRITORY)
far_nbr_list_format = FULL_LIST;
- cm_format = SYM_FULL_MATRIX;
+ cm_format = SYM_HALF_MATRIX;
#else
far_nbr_list_format = HALF_LIST;
cm_format = SYM_HALF_MATRIX;
@@ -646,8 +646,8 @@ int Init_Lists( reax_system *system, control_params *control,
//bond_top = (int*) malloc( system->total_cap * sizeof(int) );
//hb_top = (int*) malloc( system->local_cap * sizeof(int) );
- Estimate_Storages( system, control, lists,
- &Htop, hb_top, bond_top, &num_3body, comm, &matrix_dim );
+ Estimate_Storages( system, control, lists, &Htop, hb_top,
+ bond_top, &num_3body, comm, &matrix_dim, cm_format );
#if defined(NEUTRAL_TERRITORY)
Allocate_Matrix( &(workspace->H), matrix_dim, Htop, cm_format, comm );
@@ -769,8 +769,9 @@ int Init_Lists( reax_system *system, control_params *control,
bond_top = (int*) calloc( system->total_cap, sizeof(int) );
hb_top = (int*) calloc( system->local_cap, sizeof(int) );
- Estimate_Storages( system, control, lists,
- &Htop, hb_top, bond_top, &num_3body, comm, &matrix_dim );
+ //TODO: add one paramater at the end for charge matrix format - half or full
+ Estimate_Storages( system, control, lists, &Htop, hb_top,
+ bond_top, &num_3body, comm, &matrix_dim );
if ( control->hbond_cut > 0 )
{
diff --git a/PuReMD/src/linear_solvers.c b/PuReMD/src/linear_solvers.c
index f9c9345f9548a5eb031c32b9fd82bc24ea10afd3..6338c1af65fef2cdfe687d806ee8e850b4adcb6a 100644
--- a/PuReMD/src/linear_solvers.c
+++ b/PuReMD/src/linear_solvers.c
@@ -1157,9 +1157,17 @@ void Sparse_MatVec( sparse_matrix *A, real *x, real *b, int N )
si = A->start[i];
/* diagonal only contributes once */
- b[i] += A->entries[si].val * x[i];
+ if( i < A->n )
+ {
+ b[i] += A->entries[si].val * x[i];
+ k = si + 1;
+ }
+ else
+ {
+ k = si;
+ }
- for ( k = si + 1; k < A->end[i]; ++k )
+ for ( ; k < A->end[i]; ++k )
{
j = A->entries[k].j;
val = A->entries[k].val;