/*---------------------------------------------------------------------- PuReMD - Purdue ReaxFF Molecular Dynamics Program Copyright (2010) Purdue University Hasan Metin Aktulga, haktulga@cs.purdue.edu Joseph Fogarty, jcfogart@mail.usf.edu Sagar Pandit, pandit@usf.edu Ananth Y Grama, ayg@cs.purdue.edu This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details: <http://www.gnu.org/licenses/>. ----------------------------------------------------------------------*/ #if (defined(HAVE_CONFIG_H) && !defined(__CONFIG_H_)) #define __CONFIG_H_ #include "../../common/include/config.h" #endif #if defined(PURE_REAX) #include "bond_orders.h" #include "list.h" #include "vector.h" #elif defined(LAMMPS_REAX) #include "reax_bond_orders.h" #include "reax_list.h" #include "reax_vector.h" #endif #include "index_utils.h" static inline real Cf45( real p1, real p2 ) { return -EXP(-p2 / 2.0) / ( SQR( EXP(-p1 / 2.0) + EXP(p1 / 2.0) ) * (EXP(-p2 / 2.0) + EXP(p2 / 2.0)) ); } #if defined(TEST_FORCES) void Add_dBO( reax_system *system, reax_list **lists, int i, int pj, real C, rvec *v ) { int start_pj, end_pj, k; reax_list *bonds, *dBOs; bonds = lists[BONDS]; dBOs = lists[DBO]; pj = bonds->bond_list[pj].dbond_index; start_pj = Start_Index( pj, dBOs ); end_pj = End_Index( pj, dBOs ); #if defined(DEBUG_FOCUS) fprintf( stderr, "[Add_dBO] i: %d, j: %d, start: %d, end: %d\n", i, pj, start_pj, end_pj ); #endif for ( k = start_pj; k < end_pj; ++k ) { rvec_ScaledAdd( v[dBOs->dbo_list[k].wrt], C, dBOs->dbo_list[k].dBO ); } } void Add_dBOpinpi2( reax_system *system, reax_list **lists, int i, int pj, real Cpi, real Cpi2, rvec *vpi, rvec *vpi2 ) { int start_pj, end_pj, k; reax_list *bonds, *dBOs; dbond_data *dbo_k; bonds = lists[BONDS]; dBOs = lists[DBO]; pj = bonds->bond_list[pj].dbond_index; start_pj = Start_Index( pj, dBOs ); end_pj = End_Index( pj, dBOs ); #if defined(DEBUG_FOCUS) fprintf( stderr, "[Add_dBOpinpi2] i: %d, j: %d, start: %d, end: %d\n", i, pj, start_pj, end_pj ); #endif for ( k = start_pj; k < end_pj; ++k ) { dbo_k = &dBOs->dbo_list[k]; rvec_ScaledAdd( vpi[dbo_k->wrt], Cpi, dbo_k->dBOpi ); rvec_ScaledAdd( vpi2[dbo_k->wrt], Cpi2, dbo_k->dBOpi2 ); } } void Add_dBO_to_Forces( reax_system *system, reax_list **lists, int i, int pj, real C ) { reax_list *bonds; reax_list *dBOs; int start_pj, end_pj, k; bonds = lists[BONDS]; dBOs = lists[DBO]; pj = bonds->bond_list[pj].dbond_index; start_pj = Start_Index( pj, dBOs ); end_pj = End_Index( pj, dBOs ); for ( k = start_pj; k < end_pj; ++k ) { rvec_ScaledAdd( system->atoms[dBOs->dbo_list[k].wrt].f, C, dBOs->dbo_list[k].dBO ); } } void Add_dBOpinpi2_to_Forces( reax_system *system, reax_list **lists, int i, int pj, real Cpi, real Cpi2 ) { reax_list *bonds; reax_list *dBOs; dbond_data *dbo_k; int start_pj, end_pj, k; bonds = lists[BONDS]; dBOs = lists[DBO]; pj = bonds->bond_list[pj].dbond_index; start_pj = Start_Index(pj, dBOs); end_pj = End_Index(pj, dBOs); for ( k = start_pj; k < end_pj; ++k ) { dbo_k = &dBOs->dbo_list[k]; rvec_ScaledAdd( system->atoms[dbo_k->wrt].f, Cpi, dbo_k->dBOpi ); rvec_ScaledAdd( system->atoms[dbo_k->wrt].f, Cpi2, dbo_k->dBOpi2 ); } } void Add_dDelta( reax_system *system, reax_list **lists, int i, real C, rvec *v ) { int start, end, k; reax_list *dDeltas; dDeltas = lists[DDELTA]; start = Start_Index( i, dDeltas ); end = End_Index( i, dDeltas ); #if defined(DEBUG_FOCUS) fprintf( stderr, "[Add_dDelta] i: %d, start: %d, end: %d\n", i, start, end ); #endif for ( k = start; k < end; ++k ) { rvec_ScaledAdd( v[dDeltas->dDelta_list[k].wrt], C, dDeltas->dDelta_list[k].dVal ); } } void Add_dDelta_to_Forces( reax_system *system, reax_list **lists, int i, real C ) { reax_list *dDeltas; int start, end, k; dDeltas = lists[DDELTA]; start = Start_Index( i, dDeltas ); end = End_Index( i, dDeltas ); for ( k = start; k < end; ++k ) { rvec_ScaledAdd( system->atoms[dDeltas->dDelta_list[k].wrt].f, C, dDeltas->dDelta_list[k].dVal ); } } void Calculate_dBO( int i, int pj, static_storage *workspace, reax_list **lists, int *top ) { int j, k, l, start_i, end_i, end_j; reax_list *bonds, *dBOs; bond_data *nbr_l, *nbr_k; bond_order_data *bo_ij; dbond_data *top_dbo; bonds = lists[BONDS]; dBOs = lists[DBO]; j = bonds->bond_list[pj].nbr; bo_ij = &bonds->bond_list[pj].bo_data; start_i = Start_Index( i, bonds ); end_i = End_Index( i, bonds ); l = Start_Index( j, bonds ); end_j = End_Index( j, bonds ); top_dbo = &dBOs->dbo_list[ *top ]; // fprintf( stderr, "[Calculate_dBO] i: %d, pj: %d, start_i: %d, end_i: %d, start_j: %d, end_j: %d\n", i, pj, start_i, end_i, l, end_j ); for ( k = start_i; k < end_i; ++k ) { nbr_k = &bonds->bond_list[k]; for ( ; l < end_j && bonds->bond_list[l].nbr < nbr_k->nbr; ++l ) { /* These are the neighbors of j which aren't in the neighbor_list of i * Note that they might also include i! */ nbr_l = &bonds->bond_list[l]; top_dbo->wrt = nbr_l->nbr; rvec_Scale( top_dbo->dBO, -bo_ij->C3dbo, nbr_l->bo_data.dBOp ); rvec_Scale( top_dbo->dBOpi, -bo_ij->C4dbopi, nbr_l->bo_data.dBOp ); rvec_Scale( top_dbo->dBOpi2, -bo_ij->C4dbopi2, nbr_l->bo_data.dBOp ); if ( nbr_l->nbr == i ) { /* dBO */ rvec_ScaledAdd( top_dbo->dBO, bo_ij->C1dbo, bo_ij->dBOp ); rvec_ScaledAdd( top_dbo->dBO, bo_ij->C2dbo, workspace->dDeltap_self[i] ); /* dBOpi */ rvec_ScaledAdd( top_dbo->dBOpi, bo_ij->C1dbopi, bo_ij->dln_BOp_pi ); rvec_ScaledAdd( top_dbo->dBOpi, bo_ij->C2dbopi, bo_ij->dBOp ); rvec_ScaledAdd( top_dbo->dBOpi, bo_ij->C3dbopi, workspace->dDeltap_self[i] ); /* dBOpp */ rvec_ScaledAdd( top_dbo->dBOpi2, bo_ij->C1dbopi2, bo_ij->dln_BOp_pi2 ); rvec_ScaledAdd( top_dbo->dBOpi2, bo_ij->C2dbopi2, bo_ij->dBOp ); rvec_ScaledAdd( top_dbo->dBOpi2, bo_ij->C3dbopi2, workspace->dDeltap_self[i] ); } ++(*top); ++top_dbo; } /* Now we are processing neighbor k of i. */ top_dbo->wrt = nbr_k->nbr; rvec_Scale( top_dbo->dBO, -bo_ij->C2dbo, nbr_k->bo_data.dBOp ); //dBO-2 rvec_Scale( top_dbo->dBOpi, -bo_ij->C3dbopi, nbr_k->bo_data.dBOp ); //dBOpi-3 rvec_Scale( top_dbo->dBOpi2, -bo_ij->C3dbopi2, nbr_k->bo_data.dBOp );//dBOpp-3 if ( l < end_j && bonds->bond_list[l].nbr == nbr_k->nbr ) { /* common neighbor of i and j */ nbr_l = &bonds->bond_list[l]; rvec_Copy( dBOp, nbr_l->bo_data.dBOp ); rvec_ScaledAdd( top_dbo->dBO, -bo_ij->C3dbo, nbr_l->bo_data.dBOp ); //dBO,3rd rvec_ScaledAdd( top_dbo->dBOpi, -bo_ij->C4dbopi, nbr_l->bo_data.dBOp ); //dBOpi,4th rvec_ScaledAdd( top_dbo->dBOpi2, -bo_ij->C4dbopi2, nbr_l->bo_data.dBOp );//dBOpp.4th ++l; } else if ( k == pj ) { /* 1st, dBO */ rvec_ScaledAdd( top_dbo->dBO, -bo_ij->C1dbo, bo_ij->dBOp ); /* 3rd, dBO */ rvec_ScaledAdd( top_dbo->dBO, bo_ij->C3dbo, workspace->dDeltap_self[j] ); /* dBOpi, 1st */ rvec_ScaledAdd( top_dbo->dBOpi, -bo_ij->C1dbopi, bo_ij->dln_BOp_pi ); /* 2nd */ rvec_ScaledAdd( top_dbo->dBOpi, -bo_ij->C2dbopi, bo_ij->dBOp ); /* 4th */ rvec_ScaledAdd( top_dbo->dBOpi, bo_ij->C4dbopi, workspace->dDeltap_self[j] ); /* dBOpi2, 1st */ rvec_ScaledAdd( top_dbo->dBOpi2, -bo_ij->C1dbopi2, bo_ij->dln_BOp_pi2 ); /* 2nd */ rvec_ScaledAdd( top_dbo->dBOpi2, -bo_ij->C2dbopi2, bo_ij->dBOp ); /* 4th */ rvec_ScaledAdd( top_dbo->dBOpi2, bo_ij->C4dbopi2, workspace->dDeltap_self[j] ); } ++(*top); ++top_dbo; } for ( ; l < end_j; ++l ) { /* These are the remaining neighbors of j which are not in the neighbor_list of i. Note that they might also include i!*/ nbr_l = &bonds->bond_list[l]; top_dbo->wrt = nbr_l->nbr; rvec_Copy( dBOp, nbr_l->bo_data.dBOp ); rvec_Scale( top_dbo->dBO, -bo_ij->C3dbo, nbr_l->bo_data.dBOp ); //3rd, dBO rvec_Scale( top_dbo->dBOpi, -bo_ij->C4dbopi, nbr_l->bo_data.dBOp ); //4th, dBOpi rvec_Scale( top_dbo->dBOpi2, -bo_ij->C4dbopi2, nbr_l->bo_data.dBOp );//4th, dBOpp if ( nbr_l->nbr == i ) { /* dBO, 1st */ rvec_ScaledAdd( top_dbo->dBO, bo_ij->C1dbo, bo_ij->dBOp ); rvec_ScaledAdd( top_dbo->dBO, bo_ij->C2dbo, workspace->dDeltap_self[i] ); //2nd, dBO /* dBOpi, 1st */ rvec_ScaledAdd( top_dbo->dBOpi, bo_ij->C1dbopi, bo_ij->dln_BOp_pi ); rvec_ScaledAdd( top_dbo->dBOpi, bo_ij->C2dbopi, bo_ij->dBOp ); //2nd rvec_ScaledAdd( top_dbo->dBOpi, bo_ij->C3dbopi, workspace->dDeltap_self[i] ); //3rd /* dBOpipi, 1st */ rvec_ScaledAdd(top_dbo->dBOpi2, bo_ij->C1dbopi2, bo_ij->dln_BOp_pi2); rvec_ScaledAdd( top_dbo->dBOpi2, bo_ij->C2dbopi2, bo_ij->dBOp ); //2nd rvec_ScaledAdd( top_dbo->dBOpi2, bo_ij->C3dbopi2, workspace->dDeltap_self[i] );//3rd } ++(*top); ++top_dbo; } } #endif void Add_dBond_to_Forces_NPT( int i, int pj, reax_system *system, simulation_data *data, storage *workspace, reax_list **lists ) { reax_list *bonds; bond_data *nbr_j, *nbr_k; bond_order_data *bo_ij, *bo_ji; dbond_coefficients coef; rvec temp, ext_press; ivec rel_box; int pk, k, j; bonds = lists[BONDS]; nbr_j = &bonds->bond_list[pj]; j = nbr_j->nbr; bo_ij = &nbr_j->bo_data; bo_ji = &bonds->bond_list[ nbr_j->sym_index ].bo_data; coef.C1dbo = bo_ij->C1dbo * (bo_ij->Cdbo + bo_ji->Cdbo); coef.C2dbo = bo_ij->C2dbo * (bo_ij->Cdbo + bo_ji->Cdbo); coef.C3dbo = bo_ij->C3dbo * (bo_ij->Cdbo + bo_ji->Cdbo); coef.C1dbopi = bo_ij->C1dbopi * (bo_ij->Cdbopi + bo_ji->Cdbopi); coef.C2dbopi = bo_ij->C2dbopi * (bo_ij->Cdbopi + bo_ji->Cdbopi); coef.C3dbopi = bo_ij->C3dbopi * (bo_ij->Cdbopi + bo_ji->Cdbopi); coef.C4dbopi = bo_ij->C4dbopi * (bo_ij->Cdbopi + bo_ji->Cdbopi); coef.C1dbopi2 = bo_ij->C1dbopi2 * (bo_ij->Cdbopi2 + bo_ji->Cdbopi2); coef.C2dbopi2 = bo_ij->C2dbopi2 * (bo_ij->Cdbopi2 + bo_ji->Cdbopi2); coef.C3dbopi2 = bo_ij->C3dbopi2 * (bo_ij->Cdbopi2 + bo_ji->Cdbopi2); coef.C4dbopi2 = bo_ij->C4dbopi2 * (bo_ij->Cdbopi2 + bo_ji->Cdbopi2); coef.C1dDelta = bo_ij->C1dbo * (workspace->CdDelta[i] + workspace->CdDelta[j]); coef.C2dDelta = bo_ij->C2dbo * (workspace->CdDelta[i] + workspace->CdDelta[j]); coef.C3dDelta = bo_ij->C3dbo * (workspace->CdDelta[i] + workspace->CdDelta[j]); /************************************ * forces related to atom i * * first neighbors of atom i * ************************************/ for ( pk = Start_Index(i, bonds); pk < End_Index(i, bonds); ++pk ) { nbr_k = &bonds->bond_list[pk]; k = nbr_k->nbr; /* 2nd, dBO */ rvec_Scale( temp, -coef.C2dbo, nbr_k->bo_data.dBOp ); /* dDelta */ rvec_ScaledAdd( temp, -coef.C2dDelta, nbr_k->bo_data.dBOp ); /* 3rd, dBOpi */ rvec_ScaledAdd( temp, -coef.C3dbopi, nbr_k->bo_data.dBOp ); /* 3rd, dBOpi2 */ rvec_ScaledAdd( temp, -coef.C3dbopi2, nbr_k->bo_data.dBOp ); /* force */ rvec_Add( workspace->f[k], temp ); /* pressure */ rvec_iMultiply( ext_press, nbr_k->rel_box, temp ); rvec_Add( data->my_ext_press, ext_press ); } /* then atom i itself */ /* 1st, dBO */ rvec_Scale( temp, coef.C1dbo, bo_ij->dBOp ); /* 2nd, dBO */ rvec_ScaledAdd( temp, coef.C2dbo, workspace->dDeltap_self[i] ); /* 1st, dBO */ rvec_ScaledAdd( temp, coef.C1dDelta, bo_ij->dBOp ); /* 2nd, dBO */ rvec_ScaledAdd( temp, coef.C2dDelta, workspace->dDeltap_self[i] ); /* 1st, dBOpi */ rvec_ScaledAdd( temp, coef.C1dbopi, bo_ij->dln_BOp_pi ); /* 2nd, dBOpi */ rvec_ScaledAdd( temp, coef.C2dbopi, bo_ij->dBOp ); /* 3rd, dBOpi */ rvec_ScaledAdd( temp, coef.C3dbopi, workspace->dDeltap_self[i] ); /* 1st, dBO_pi2 */ rvec_ScaledAdd( temp, coef.C1dbopi2, bo_ij->dln_BOp_pi2 ); /* 2nd, dBO_pi2 */ rvec_ScaledAdd( temp, coef.C2dbopi2, bo_ij->dBOp ); /* 3rd, dBO_pi2 */ rvec_ScaledAdd( temp, coef.C3dbopi2, workspace->dDeltap_self[i] ); /* force */ rvec_Add( workspace->f[i], temp ); /* ext pressure due to i is dropped, counting force on j will be enough */ /****************************************************** * forces and pressure related to atom j * * first neighbors of atom j * ******************************************************/ for ( pk = Start_Index(j, bonds); pk < End_Index(j, bonds); ++pk ) { nbr_k = &bonds->bond_list[pk]; k = nbr_k->nbr; /* 3rd, dBO */ rvec_Scale( temp, -coef.C3dbo, nbr_k->bo_data.dBOp ); /* dDelta */ rvec_ScaledAdd( temp, -coef.C3dDelta, nbr_k->bo_data.dBOp ); /* 4th, dBOpi */ rvec_ScaledAdd( temp, -coef.C4dbopi, nbr_k->bo_data.dBOp ); /* 4th, dBOpi2 */ rvec_ScaledAdd( temp, -coef.C4dbopi2, nbr_k->bo_data.dBOp ); /* force */ rvec_Add( workspace->f[k], temp ); /* pressure */ if ( k != i ) { ivec_Sum( rel_box, nbr_k->rel_box, nbr_j->rel_box ); //rel_box(k, i) rvec_iMultiply( ext_press, rel_box, temp ); rvec_Add( data->my_ext_press, ext_press ); } } /* then atom j itself */ /* 1st, dBO */ rvec_Scale( temp, -coef.C1dbo, bo_ij->dBOp ); /* 2nd, dBO */ rvec_ScaledAdd( temp, coef.C3dbo, workspace->dDeltap_self[j] ); /* 1st, dBO */ rvec_ScaledAdd( temp, -coef.C1dDelta, bo_ij->dBOp ); /* 2nd, dBO */ rvec_ScaledAdd( temp, coef.C3dDelta, workspace->dDeltap_self[j] ); /* 1st, dBOpi */ rvec_ScaledAdd( temp, -coef.C1dbopi, bo_ij->dln_BOp_pi ); /* 2nd, dBOpi */ rvec_ScaledAdd( temp, -coef.C2dbopi, bo_ij->dBOp ); /* 3rd, dBOpi */ rvec_ScaledAdd( temp, coef.C4dbopi, workspace->dDeltap_self[j] ); /* 1st, dBOpi2 */ rvec_ScaledAdd( temp, -coef.C1dbopi2, bo_ij->dln_BOp_pi2 ); /* 2nd, dBOpi2 */ rvec_ScaledAdd( temp, -coef.C2dbopi2, bo_ij->dBOp ); /* 3rd, dBOpi2 */ rvec_ScaledAdd( temp, coef.C4dbopi2, workspace->dDeltap_self[j] ); /* force */ rvec_Add( workspace->f[j], temp ); /* pressure */ rvec_iMultiply( ext_press, nbr_j->rel_box, temp ); rvec_Add( data->my_ext_press, ext_press ); } void Add_dBond_to_Forces( int i, int pj, reax_system *system, simulation_data *data, storage *workspace, reax_list **lists ) { reax_list *bonds; bond_data *nbr_j, *nbr_k; bond_order_data *bo_ij, *bo_ji; dbond_coefficients coef; int pk, k, j; bonds = lists[BONDS]; nbr_j = &bonds->bond_list[pj]; j = nbr_j->nbr; bo_ij = &nbr_j->bo_data; bo_ji = &bonds->bond_list[ nbr_j->sym_index ].bo_data; coef.C1dbo = bo_ij->C1dbo * (bo_ij->Cdbo + bo_ji->Cdbo); coef.C2dbo = bo_ij->C2dbo * (bo_ij->Cdbo + bo_ji->Cdbo); coef.C3dbo = bo_ij->C3dbo * (bo_ij->Cdbo + bo_ji->Cdbo); coef.C1dbopi = bo_ij->C1dbopi * (bo_ij->Cdbopi + bo_ji->Cdbopi); coef.C2dbopi = bo_ij->C2dbopi * (bo_ij->Cdbopi + bo_ji->Cdbopi); coef.C3dbopi = bo_ij->C3dbopi * (bo_ij->Cdbopi + bo_ji->Cdbopi); coef.C4dbopi = bo_ij->C4dbopi * (bo_ij->Cdbopi + bo_ji->Cdbopi); coef.C1dbopi2 = bo_ij->C1dbopi2 * (bo_ij->Cdbopi2 + bo_ji->Cdbopi2); coef.C2dbopi2 = bo_ij->C2dbopi2 * (bo_ij->Cdbopi2 + bo_ji->Cdbopi2); coef.C3dbopi2 = bo_ij->C3dbopi2 * (bo_ij->Cdbopi2 + bo_ji->Cdbopi2); coef.C4dbopi2 = bo_ij->C4dbopi2 * (bo_ij->Cdbopi2 + bo_ji->Cdbopi2); coef.C1dDelta = bo_ij->C1dbo * (workspace->CdDelta[i] + workspace->CdDelta[j]); coef.C2dDelta = bo_ij->C2dbo * (workspace->CdDelta[i] + workspace->CdDelta[j]); coef.C3dDelta = bo_ij->C3dbo * (workspace->CdDelta[i] + workspace->CdDelta[j]); for ( pk = Start_Index(i, bonds); pk < End_Index(i, bonds); ++pk ) { nbr_k = &bonds->bond_list[pk]; k = nbr_k->nbr; /* 2nd, dBO */ rvec_ScaledAdd( workspace->f[k], -coef.C2dbo, nbr_k->bo_data.dBOp ); /* dDelta */ rvec_ScaledAdd( workspace->f[k], -coef.C2dDelta, nbr_k->bo_data.dBOp ); /* 3rd, dBOpi */ rvec_ScaledAdd( workspace->f[k], -coef.C3dbopi, nbr_k->bo_data.dBOp ); /* 3rd, dBOpi2 */ rvec_ScaledAdd( workspace->f[k], -coef.C3dbopi2, nbr_k->bo_data.dBOp ); } /* 1st, dBO */ rvec_ScaledAdd( workspace->f[i], coef.C1dbo, bo_ij->dBOp ); /* 2nd, dBO */ rvec_ScaledAdd( workspace->f[i], coef.C2dbo, workspace->dDeltap_self[i] ); /* 1st, dBO */ rvec_ScaledAdd( workspace->f[i], coef.C1dDelta, bo_ij->dBOp ); /* 2nd, dBO */ rvec_ScaledAdd( workspace->f[i], coef.C2dDelta, workspace->dDeltap_self[i] ); /* 1st, dBOpi */ rvec_ScaledAdd( workspace->f[i], coef.C1dbopi, bo_ij->dln_BOp_pi ); /* 2nd, dBOpi */ rvec_ScaledAdd( workspace->f[i], coef.C2dbopi, bo_ij->dBOp ); /* 3rd, dBOpi */ rvec_ScaledAdd( workspace->f[i], coef.C3dbopi, workspace->dDeltap_self[i] ); /* 1st, dBO_pi2 */ rvec_ScaledAdd( workspace->f[i], coef.C1dbopi2, bo_ij->dln_BOp_pi2 ); /* 2nd, dBO_pi2 */ rvec_ScaledAdd( workspace->f[i], coef.C2dbopi2, bo_ij->dBOp ); /* 3rd, dBO_pi2 */ rvec_ScaledAdd( workspace->f[i], coef.C3dbopi2, workspace->dDeltap_self[i] ); for ( pk = Start_Index(j, bonds); pk < End_Index(j, bonds); ++pk ) { nbr_k = &bonds->bond_list[pk]; k = nbr_k->nbr; /* 3rd, dBO */ rvec_ScaledAdd( workspace->f[k], -coef.C3dbo, nbr_k->bo_data.dBOp ); /* dDelta */ rvec_ScaledAdd( workspace->f[k], -coef.C3dDelta, nbr_k->bo_data.dBOp ); /* 4th, dBOpi */ rvec_ScaledAdd( workspace->f[k], -coef.C4dbopi, nbr_k->bo_data.dBOp ); /* 4th, dBOpi2 */ rvec_ScaledAdd( workspace->f[k], -coef.C4dbopi2, nbr_k->bo_data.dBOp ); } /* 1st, dBO */ rvec_ScaledAdd( workspace->f[j], -coef.C1dbo, bo_ij->dBOp ); /* 2nd, dBO */ rvec_ScaledAdd( workspace->f[j], coef.C3dbo, workspace->dDeltap_self[j] ); /* 1st, dBO */ rvec_ScaledAdd( workspace->f[j], -coef.C1dDelta, bo_ij->dBOp ); /* 2nd, dBO */ rvec_ScaledAdd( workspace->f[j], coef.C3dDelta, workspace->dDeltap_self[j] ); /* 1st, dBOpi */ rvec_ScaledAdd( workspace->f[j], -coef.C1dbopi, bo_ij->dln_BOp_pi ); /* 2nd, dBOpi */ rvec_ScaledAdd( workspace->f[j], -coef.C2dbopi, bo_ij->dBOp ); /* 3rd, dBOpi */ rvec_ScaledAdd( workspace->f[j], coef.C4dbopi, workspace->dDeltap_self[j] ); /* 1st, dBOpi2 */ rvec_ScaledAdd( workspace->f[j], -coef.C1dbopi2, bo_ij->dln_BOp_pi2 ); /* 2nd, dBOpi2 */ rvec_ScaledAdd( workspace->f[j], -coef.C2dbopi2, bo_ij->dBOp ); /* 3rd, dBOpi2 */ rvec_ScaledAdd( workspace->f[j], coef.C4dbopi2, workspace->dDeltap_self[j] ); } static inline void Copy_Bond_Order_Data( bond_order_data *dest, bond_order_data *src ) { dest->BO = src->BO; dest->BO_s = src->BO_s; dest->BO_pi = src->BO_pi; dest->BO_pi2 = src->BO_pi2; rvec_Scale( dest->dBOp, -1.0, src->dBOp ); rvec_Scale( dest->dln_BOp_s, -1.0, src->dln_BOp_s ); rvec_Scale( dest->dln_BOp_pi, -1.0, src->dln_BOp_pi ); rvec_Scale( dest->dln_BOp_pi2, -1.0, src->dln_BOp_pi2 ); } /* Compute the bond order term between atoms i and j, * and if this term exceeds the cutoff bo_cut, then adds * BOTH atoms the bonds list (i.e., compute term once * and copy to avoid redundant computation) */ int BOp( storage * const workspace, reax_list * const bonds, real bo_cut, int i, int btop_i, int j, ivec * const rel_box, real d, rvec * const dvec, int far_nbr_list_format, single_body_parameters const * const sbp_i, single_body_parameters const * const sbp_j, two_body_parameters const * const twbp ) { real r2, C12, C34, C56; real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2; real BO, BO_s, BO_pi, BO_pi2; bond_data *ibond, *jbond; bond_order_data *bo_ij, *bo_ji; int btop_j; r2 = SQR( d ); if ( sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0 ) { C12 = twbp->p_bo1 * POW( d / twbp->r_s, twbp->p_bo2 ); BO_s = (1.0 + bo_cut) * exp( C12 ); } else { C12 = 0.0; BO_s = 0.0; } if ( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0 ) { C34 = twbp->p_bo3 * POW( d / twbp->r_p, twbp->p_bo4 ); BO_pi = exp( C34 ); } else { C34 = 0.0; BO_pi = 0.0; } if ( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0 ) { C56 = twbp->p_bo5 * POW( d / twbp->r_pp, twbp->p_bo6 ); BO_pi2 = exp( C56 ); } else { C56 = 0.0; BO_pi2 = 0.0; } /* Initially BO values are the uncorrected ones, page 1 */ BO = BO_s + BO_pi + BO_pi2; if ( BO >= bo_cut ) { /****** bonds i-j and j-i ******/ ibond = &bonds->bond_list[btop_i]; btop_j = End_Index( j, bonds ); jbond = &bonds->bond_list[btop_j]; ibond->nbr = j; ibond->d = d; rvec_Copy( ibond->dvec, *dvec ); ivec_Copy( ibond->rel_box, *rel_box ); ibond->dbond_index = btop_i; ibond->sym_index = btop_j; jbond->nbr = i; jbond->d = d; rvec_Scale( jbond->dvec, -1.0, *dvec ); ivec_Scale( jbond->rel_box, -1.0, *rel_box ); jbond->dbond_index = btop_i; jbond->sym_index = btop_i; Set_End_Index( j, btop_j + 1, bonds ); bo_ij = &ibond->bo_data; bo_ij->BO = BO; bo_ij->BO_s = BO_s; bo_ij->BO_pi = BO_pi; bo_ij->BO_pi2 = BO_pi2; bo_ji = &jbond->bo_data; bo_ji->BO = BO; bo_ji->BO_s = BO_s; bo_ji->BO_pi = BO_pi; bo_ji->BO_pi2 = BO_pi2; /* Bond Order page2-3, derivative of total bond order prime */ Cln_BOp_s = twbp->p_bo2 * C12 / r2; Cln_BOp_pi = twbp->p_bo4 * C34 / r2; Cln_BOp_pi2 = twbp->p_bo6 * C56 / r2; /* Only dln_BOp_xx wrt. dr_i is stored here, note that * dln_BOp_xx/dr_i = -dln_BOp_xx/dr_j and all others are 0 */ rvec_Scale( bo_ij->dln_BOp_s, -1.0 * bo_ij->BO_s * Cln_BOp_s, ibond->dvec ); rvec_Scale( bo_ij->dln_BOp_pi, -1.0 * bo_ij->BO_pi * Cln_BOp_pi, ibond->dvec ); rvec_Scale( bo_ij->dln_BOp_pi2, -1.0 * bo_ij->BO_pi2 * Cln_BOp_pi2, ibond->dvec ); rvec_Scale( bo_ji->dln_BOp_s, -1.0, bo_ij->dln_BOp_s ); rvec_Scale( bo_ji->dln_BOp_pi, -1.0, bo_ij->dln_BOp_pi ); rvec_Scale( bo_ji->dln_BOp_pi2, -1.0, bo_ij->dln_BOp_pi2 ); /* Only dBOp wrt. dr_i is stored here, note that * dBOp/dr_i = -dBOp/dr_j and all others are 0 */ rvec_Scale( bo_ij->dBOp, -1.0 * (bo_ij->BO_s * Cln_BOp_s + bo_ij->BO_pi * Cln_BOp_pi + bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec ); rvec_Scale( bo_ji->dBOp, -1.0, bo_ij->dBOp ); rvec_Add( workspace->dDeltap_self[i], bo_ij->dBOp ); rvec_Add( workspace->dDeltap_self[j], bo_ji->dBOp ); bo_ij->BO_s -= bo_cut; bo_ij->BO -= bo_cut; /* currently total_BOp */ workspace->total_bond_order[i] += bo_ij->BO; bo_ij->Cdbo = 0.0; bo_ij->Cdbopi = 0.0; bo_ij->Cdbopi2 = 0.0; bo_ji->BO_s -= bo_cut; bo_ji->BO -= bo_cut; workspace->total_bond_order[j] += bo_ji->BO; //currently total_BOp bo_ji->Cdbo = 0.0; bo_ji->Cdbopi = 0.0; bo_ji->Cdbopi2 = 0.0; return TRUE; } return FALSE; } /* Compute the bond order term between atoms i and j, * and if this term exceeds the cutoff bo_cut, then adds * to the bond list according to the following convention: * * if the far neighbor list is stored in half format, * add BOTH atoms to each other's portion of the bond list * * if the far neighbor list is stored in full format, * add atom i to atom j's bonds list ONLY */ int BOp_redundant( storage *workspace, reax_list *bonds, real bo_cut, int i, int btop_i, int j, ivec *rel_box, real d, rvec *dvec, int far_nbr_list_format, single_body_parameters *sbp_i, single_body_parameters *sbp_j, two_body_parameters *twbp ) { real r2, C12, C34, C56; real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2; real BO, BO_s, BO_pi, BO_pi2; bond_data *ibond; bond_order_data *bo_ij; int btop_j; bond_data *jbond; bond_order_data *bo_ji; r2 = SQR(d); if ( sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0 ) { C12 = twbp->p_bo1 * pow( d / twbp->r_s, twbp->p_bo2 ); BO_s = (1.0 + bo_cut) * exp( C12 ); } else { C12 = 0.0; BO_s = 0.0; } if ( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0 ) { C34 = twbp->p_bo3 * pow( d / twbp->r_p, twbp->p_bo4 ); BO_pi = exp( C34 ); } else { C34 = 0.0; BO_pi = 0.0; } if ( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0 ) { C56 = twbp->p_bo5 * pow( d / twbp->r_pp, twbp->p_bo6 ); BO_pi2 = exp( C56 ); } else { C56 = 0.0; BO_pi2 = 0.0; } /* Initially BO values are the uncorrected ones, page 1 */ BO = BO_s + BO_pi + BO_pi2; if ( BO >= bo_cut ) { /****** bonds i-j and j-i ******/ ibond = &bonds->bond_list[btop_i]; if ( far_nbr_list_format == HALF_LIST ) { btop_j = End_Index( j, bonds ); jbond = &bonds->bond_list[btop_j]; } ibond->nbr = j; ibond->d = d; rvec_Copy( ibond->dvec, *dvec ); ivec_Copy( ibond->rel_box, *rel_box ); ibond->dbond_index = btop_i; if ( far_nbr_list_format == HALF_LIST ) { ibond->sym_index = btop_j; jbond->nbr = i; jbond->d = d; rvec_Scale( jbond->dvec, -1.0, *dvec ); ivec_Scale( jbond->rel_box, -1.0, *rel_box ); jbond->dbond_index = btop_i; jbond->sym_index = btop_i; Set_End_Index( j, btop_j + 1, bonds ); } bo_ij = &ibond->bo_data; bo_ij->BO = BO; bo_ij->BO_s = BO_s; bo_ij->BO_pi = BO_pi; bo_ij->BO_pi2 = BO_pi2; if ( far_nbr_list_format == HALF_LIST ) { bo_ji = &jbond->bo_data; bo_ji->BO = BO; bo_ji->BO_s = BO_s; bo_ji->BO_pi = BO_pi; bo_ji->BO_pi2 = BO_pi2; } /* Bond Order page2-3, derivative of total bond order prime */ Cln_BOp_s = twbp->p_bo2 * C12 / r2; Cln_BOp_pi = twbp->p_bo4 * C34 / r2; Cln_BOp_pi2 = twbp->p_bo6 * C56 / r2; /* Only dln_BOp_xx wrt. dr_i is stored here, note that * dln_BOp_xx/dr_i = -dln_BOp_xx/dr_j and all others are 0 */ rvec_Scale( bo_ij->dln_BOp_s, -1.0 * bo_ij->BO_s * Cln_BOp_s, ibond->dvec ); rvec_Scale( bo_ij->dln_BOp_pi, -1.0 * bo_ij->BO_pi * Cln_BOp_pi, ibond->dvec ); rvec_Scale( bo_ij->dln_BOp_pi2, -1.0 * bo_ij->BO_pi2 * Cln_BOp_pi2, ibond->dvec ); if ( far_nbr_list_format == HALF_LIST ) { rvec_Scale( bo_ji->dln_BOp_s, -1.0, bo_ij->dln_BOp_s ); rvec_Scale( bo_ji->dln_BOp_pi, -1.0, bo_ij->dln_BOp_pi ); rvec_Scale( bo_ji->dln_BOp_pi2, -1.0, bo_ij->dln_BOp_pi2 ); } /* Only dBOp wrt. dr_i is stored here, note that * dBOp/dr_i = -dBOp/dr_j and all others are 0 */ rvec_Scale( bo_ij->dBOp, -1.0 * (bo_ij->BO_s * Cln_BOp_s + bo_ij->BO_pi * Cln_BOp_pi + bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec ); if ( far_nbr_list_format == HALF_LIST ) { rvec_Scale( bo_ji->dBOp, -1.0, bo_ij->dBOp ); } rvec_Add( workspace->dDeltap_self[i], bo_ij->dBOp ); if ( far_nbr_list_format == HALF_LIST ) { rvec_Add( workspace->dDeltap_self[j], bo_ji->dBOp ); } bo_ij->BO_s -= bo_cut; bo_ij->BO -= bo_cut; workspace->total_bond_order[i] += bo_ij->BO; //currently total_BOp bo_ij->Cdbo = 0.0; bo_ij->Cdbopi = 0.0; bo_ij->Cdbopi2 = 0.0; if ( far_nbr_list_format == HALF_LIST ) { bo_ji->BO_s -= bo_cut; bo_ji->BO -= bo_cut; workspace->total_bond_order[j] += bo_ji->BO; //currently total_BOp bo_ji->Cdbo = 0.0; bo_ji->Cdbopi = 0.0; bo_ji->Cdbopi2 = 0.0; } return TRUE; } return FALSE; } /* A very important and crucial assumption here is that each segment * belonging to a different atom in nbrhoods->nbr_list is sorted in its own. * This can either be done in the general coordinator function or here */ void BO( reax_system *system, control_params *control, simulation_data *data, storage *workspace, reax_list **lists, output_controls *out_control ) { int i, j, pj, type_i, type_j; int start_i, end_i; int sym_index; real val_i, Deltap_i, Deltap_boc_i; real val_j, Deltap_j, Deltap_boc_j; real f1, f2, f3, f4, f5, f4f5, exp_f4, exp_f5; real exp_p1i, exp_p2i, exp_p1j, exp_p2j; real temp, u1_ij, u1_ji, Cf1A_ij, Cf1B_ij, Cf1_ij, Cf1_ji; real Cf45_ij, Cf45_ji, p_lp1; real A0_ij, A1_ij, A2_ij, A2_ji, A3_ij, A3_ji; real explp1, p_boc1, p_boc2; single_body_parameters *sbp_i, *sbp_j; two_body_parameters *twbp; bond_order_data *bo_ij, *bo_ji; reax_list *bond_list; #if defined(TEST_FORCES) int k, pk, start_j, end_j; int top_dbo, top_dDelta; dbond_data *pdbo; dDelta_data *ptop_dDelta; reax_list *dDeltas, *dBOs; #endif p_lp1 = system->reax_param.gp.l[15]; p_boc1 = system->reax_param.gp.l[0]; p_boc2 = system->reax_param.gp.l[1]; bond_list = lists[BONDS]; #if defined(TEST_FORCES) top_dbo = 0; top_dDelta = 0; dDeltas = lists[DDELTAS]; dBOs = lists[DBOS]; #endif /* Calculate Deltaprime, Deltaprime_boc values */ for ( i = 0; i < system->N; ++i ) { type_i = system->my_atoms[i].type; sbp_i = &system->reax_param.sbp[type_i]; workspace->Deltap[i] = workspace->total_bond_order[i] - sbp_i->valency; workspace->Deltap_boc[i] = workspace->total_bond_order[i] - sbp_i->valency_val; workspace->total_bond_order[i] = 0.0; } /* Corrected Bond Order calculations */ for ( i = 0; i < system->N; ++i ) { type_i = system->my_atoms[i].type; sbp_i = &system->reax_param.sbp[type_i]; val_i = sbp_i->valency; Deltap_i = workspace->Deltap[i]; Deltap_boc_i = workspace->Deltap_boc[i]; start_i = Start_Index( i, bond_list ); end_i = End_Index( i, bond_list ); for ( pj = start_i; pj < end_i; ++pj ) { j = bond_list->bond_list[pj].nbr; type_j = system->my_atoms[j].type; bo_ij = &bond_list->bond_list[pj].bo_data; if ( i < j || workspace->bond_mark[j] > 3 ) { twbp = &system->reax_param.tbp[ index_tbp(type_i, type_j, system->reax_param.num_atom_types) ]; #if defined(TEST_FORCES) Set_Start_Index( pj, top_dbo, dBOs ); #endif if ( twbp->ovc < 0.001 && twbp->v13cor < 0.001 ) { /* There is no correction to bond orders nor to derivatives of * bond order prime! So we leave bond orders unchanged and * set derivative of bond order coefficients s.t. * dBO = dBOp & dBOxx = dBOxxp in Add_dBO_to_Forces */ bo_ij->C1dbo = 1.0; bo_ij->C2dbo = 0.0; bo_ij->C3dbo = 0.0; bo_ij->C1dbopi = 1.0; bo_ij->C2dbopi = 0.0; bo_ij->C3dbopi = 0.0; bo_ij->C4dbopi = 0.0; bo_ij->C1dbopi2 = 1.0; bo_ij->C2dbopi2 = 0.0; bo_ij->C3dbopi2 = 0.0; bo_ij->C4dbopi2 = 0.0; #if defined(TEST_FORCES) pdbo = &dBOs->dbo_list[ top_dbo ]; /* compute dBO_ij/dr_i */ pdbo->wrt = i; rvec_Copy( pdbo->dBO, bo_ij->dBOp ); rvec_Scale( pdbo->dBOpi, bo_ij->BO_pi, bo_ij->dln_BOp_pi ); rvec_Scale( pdbo->dBOpi2, bo_ij->BO_pi2, bo_ij->dln_BOp_pi2); /* compute dBO_ij/dr_j */ pdbo = &dBOs->dbo_list[ top_dbo + 1 ]; pdbo->wrt = j; rvec_Scale( pdbo->dBO, -1.0, bo_ij->dBOp ); rvec_Scale( pdbo->dBOpi, -bo_ij->BO_pi, bo_ij->dln_BOp_pi ); rvec_Scale(pdbo->dBOpi2, -bo_ij->BO_pi2, bo_ij->dln_BOp_pi2); top_dbo += 2; #endif } else { val_j = system->reax_param.sbp[type_j].valency; Deltap_j = workspace->Deltap[j]; Deltap_boc_j = workspace->Deltap_boc[j]; /* on page 1 */ if ( twbp->ovc >= 0.001 ) { /* Correction for overcoordination */ exp_p1i = EXP( -p_boc1 * Deltap_i ); exp_p2i = EXP( -p_boc2 * Deltap_i ); exp_p1j = EXP( -p_boc1 * Deltap_j ); exp_p2j = EXP( -p_boc2 * Deltap_j ); f2 = exp_p1i + exp_p1j; f3 = -1.0 / p_boc2 * LOG( 0.5 * ( exp_p2i + exp_p2j ) ); f1 = 0.5 * ( ( val_i + f2 ) / ( val_i + f2 + f3 ) + ( val_j + f2 ) / ( val_j + f2 + f3 ) ); /* Now come the derivates */ /* Bond Order pages 5-7, derivative of f1 */ temp = f2 + f3; u1_ij = val_i + temp; u1_ji = val_j + temp; Cf1A_ij = 0.5 * f3 * (1.0 / SQR( u1_ij ) + 1.0 / SQR( u1_ji )); Cf1B_ij = -0.5 * (( u1_ij - f3 ) / SQR( u1_ij ) + ( u1_ji - f3 ) / SQR( u1_ji )); //Cf1_ij = -Cf1A_ij * p_boc1 * exp_p1i + // Cf1B_ij * exp_p2i / ( exp_p2i + exp_p2j ); Cf1_ij = 0.50 * ( -p_boc1 * exp_p1i / u1_ij - ((val_i + f2) / SQR(u1_ij)) * ( -p_boc1 * exp_p1i + exp_p2i / ( exp_p2i + exp_p2j ) ) + -p_boc1 * exp_p1i / u1_ji - ((val_j + f2) / SQR(u1_ji)) * ( -p_boc1 * exp_p1i + exp_p2i / ( exp_p2i + exp_p2j ) )); Cf1_ji = -Cf1A_ij * p_boc1 * exp_p1j + Cf1B_ij * exp_p2j / ( exp_p2i + exp_p2j ); } else { /* No overcoordination correction! */ f1 = 1.0; Cf1_ij = 0.0; Cf1_ji = 0.0; } if ( twbp->v13cor >= 0.001 ) { /* Correction for 1-3 bond orders */ exp_f4 = EXP( -twbp->p_boc3 * (twbp->p_boc4 * SQR( bo_ij->BO ) - Deltap_boc_i) + twbp->p_boc5 ); exp_f5 = EXP( -twbp->p_boc3 * (twbp->p_boc4 * SQR( bo_ij->BO ) - Deltap_boc_j) + twbp->p_boc5 ); f4 = 1.0 / (1.0 + exp_f4); f5 = 1.0 / (1.0 + exp_f5); f4f5 = f4 * f5; /* Bond Order pages 8-9, derivative of f4 and f5 */ // temp = twbp->p_boc5 // - twbp->p_boc3 * twbp->p_boc4 * SQR( bo_ij->BO ); // u_ij = temp + twbp->p_boc3 * Deltap_boc_i; // u_ji = temp + twbp->p_boc3 * Deltap_boc_j; // Cf45_ij = Cf45( u_ij, u_ji ) / f4f5; // Cf45_ji = Cf45( u_ji, u_ij ) / f4f5; Cf45_ij = -f4 * exp_f4; Cf45_ji = -f5 * exp_f5; } else { f4 = 1.0; f5 = 1.0; f4f5 = 1.0; Cf45_ij = 0.0; Cf45_ji = 0.0; } /* Bond Order page 10, derivative of total bond order */ A0_ij = f1 * f4f5; A1_ij = -2.0 * twbp->p_boc3 * twbp->p_boc4 * bo_ij->BO * (Cf45_ij + Cf45_ji); A2_ij = Cf1_ij / f1 + twbp->p_boc3 * Cf45_ij; A2_ji = Cf1_ji / f1 + twbp->p_boc3 * Cf45_ji; A3_ij = A2_ij + Cf1_ij / f1; A3_ji = A2_ji + Cf1_ji / f1; /* find corrected bond order values and their deriv coefs */ bo_ij->BO = bo_ij->BO * A0_ij; bo_ij->BO_pi = bo_ij->BO_pi * A0_ij * f1; bo_ij->BO_pi2 = bo_ij->BO_pi2 * A0_ij * f1; bo_ij->BO_s = bo_ij->BO - ( bo_ij->BO_pi + bo_ij->BO_pi2 ); bo_ij->C1dbo = A0_ij + bo_ij->BO * A1_ij; bo_ij->C2dbo = bo_ij->BO * A2_ij; bo_ij->C3dbo = bo_ij->BO * A2_ji; bo_ij->C1dbopi = f1 * f1 * f4 * f5; bo_ij->C2dbopi = bo_ij->BO_pi * A1_ij; bo_ij->C3dbopi = bo_ij->BO_pi * A3_ij; bo_ij->C4dbopi = bo_ij->BO_pi * A3_ji; bo_ij->C1dbopi2 = f1 * f1 * f4 * f5; bo_ij->C2dbopi2 = bo_ij->BO_pi2 * A1_ij; bo_ij->C3dbopi2 = bo_ij->BO_pi2 * A3_ij; bo_ij->C4dbopi2 = bo_ij->BO_pi2 * A3_ji; #if defined(TEST_FORCES) Calculate_dBO( i, pj, workspace, lists, &top_dbo ); #endif } /* neglect weak bonds */ if ( bo_ij->BO < 1.0e-10 ) { bo_ij->BO = 0.0; } if ( bo_ij->BO_s < 1.0e-10 ) { bo_ij->BO_s = 0.0; } if ( bo_ij->BO_pi < 1.0e-10 ) { bo_ij->BO_pi = 0.0; } if ( bo_ij->BO_pi2 < 1.0e-10 ) { bo_ij->BO_pi2 = 0.0; } /* now keeps total_BO */ workspace->total_bond_order[i] += bo_ij->BO; #if defined(TEST_FORCES) Set_End_Index( pj, top_dbo, dBOs ); Add_dBO( system, lists, i, pj, 1.0, workspace->dDelta ); #endif } else { /* We only need to update bond orders from bo_ji * everything else is set in uncorrected_bo calculations */ sym_index = bond_list->bond_list[pj].sym_index; bo_ji = &bond_list->bond_list[ sym_index ].bo_data; bo_ij->BO = bo_ji->BO; bo_ij->BO_s = bo_ji->BO_s; bo_ij->BO_pi = bo_ji->BO_pi; bo_ij->BO_pi2 = bo_ji->BO_pi2; /* now keeps total_BO */ workspace->total_bond_order[i] += bo_ij->BO; #if defined(TEST_FORCES) Add_dBO( system, lists, j, sym_index, 1.0, workspace->dDelta ); #endif } } #if defined(TEST_FORCES) Set_Start_Index( i, top_dDelta, dDeltas ); ptop_dDelta = &dDeltas->dDelta_list[top_dDelta]; for ( pj = start_i; pj < end_i; ++pj ) { j = bond_list->bond_list[pj].nbr; if ( !rvec_isZero( workspace->dDelta[j] ) ) { ptop_dDelta->wrt = j; rvec_Copy( ptop_dDelta->dVal, workspace->dDelta[j] ); rvec_MakeZero( workspace->dDelta[j] ); ++top_dDelta; ++ptop_dDelta; } start_j = Start_Index( j, bond_list ); end_j = End_Index( j, bond_list ); for ( pk = start_j; pk < end_j; ++pk ) { k = bond_list->bond_list[pk].nbr; if ( !rvec_isZero( workspace->dDelta[k] ) ) { ptop_dDelta->wrt = k; rvec_Copy( ptop_dDelta->dVal, workspace->dDelta[k] ); rvec_MakeZero( workspace->dDelta[k] ); ++top_dDelta; ++ptop_dDelta; } } } Set_End_Index( i, top_dDelta, dDeltas ); #endif } /* Calculate some helper variables that are used at many places * throughout force calculations */ for ( i = 0; i < system->N; ++i ) { type_j = system->my_atoms[i].type; sbp_j = &system->reax_param.sbp[ type_j ]; workspace->Delta[i] = workspace->total_bond_order[i] - sbp_j->valency; workspace->Delta_e[i] = workspace->total_bond_order[i] - sbp_j->valency_e; workspace->Delta_boc[i] = workspace->total_bond_order[i] - sbp_j->valency_boc; workspace->vlpex[i] = workspace->Delta_e[i] - 2.0 * (int)(workspace->Delta_e[i] / 2.0); explp1 = EXP(-p_lp1 * SQR(2.0 + workspace->vlpex[i])); workspace->nlp[i] = explp1 - (int)(workspace->Delta_e[i] / 2.0); workspace->Delta_lp[i] = sbp_j->nlp_opt - workspace->nlp[i]; workspace->Clp[i] = 2.0 * p_lp1 * explp1 * (2.0 + workspace->vlpex[i]); /* Adri uses different dDelta_lp values than the ones in notes... */ workspace->dDelta_lp[i] = workspace->Clp[i]; // workspace->dDelta_lp[i] = workspace->Clp[i] + (0.5 - workspace->Clp[i]) // * ((FABS(workspace->Delta_e[i] / 2.0 // - (int)(workspace->Delta_e[i] / 2.0)) < 0.1) ? 1 : 0 ); if ( sbp_j->mass > 21.0 ) { workspace->nlp_temp[i] = 0.5 * (sbp_j->valency_e - sbp_j->valency); workspace->Delta_lp_temp[i] = sbp_j->nlp_opt - workspace->nlp_temp[i]; workspace->dDelta_lp_temp[i] = 0.0; } else { workspace->nlp_temp[i] = workspace->nlp[i]; workspace->Delta_lp_temp[i] = sbp_j->nlp_opt - workspace->nlp_temp[i]; workspace->dDelta_lp_temp[i] = workspace->Clp[i]; } } #if defined(TEST_ENERGIES) || defined(TEST_FORCES) Print_Bond_List( system, control, data, lists, out_control); #endif }