/* */ /* CDDL HEADER START */ /* */ /* The contents of this file are subject to the terms of the Common */ /* Development and Distribution License Version 1.0 (the "License"). */ /* */ /* You can obtain a copy of the license at */ /* http://www.opensource.org/licenses/CDDL-1.0. See the License for the */ /* specific language governing permissions and limitations under the License. */ /* */ /* When distributing Covered Code, include this CDDL HEADER in each file and */ /* include the License file in a prominent location with the name */ /* LICENSE.CDDL. If applicable, add the following below this CDDL HEADER, */ /* with the fields enclosed by brackets "[]" replaced with your own */ /* identifying information: */ /* */ /* Portions Copyright (c) [yyyy] [name of copyright owner]. */ /* All rights reserved. */ /* */ /* CDDL HEADER END */ /* */ /* */ /* Copyright (c) 2019, Regents of the University of Minnesota. */ /* All rights reserved. */ /* */ /* Contributors: */ /* Daniel S. Karls */ /* Ellad B. Tadmor */ /* Ryan S. Elliott */ /* */ #include "KIM_LogMacros.h" #include "KIM_ModelDriverHeaders.h" #include #include #include #include #include "cluster.inc" #define TRUE 1 #define FALSE 0 /******************************************************************************/ /* Below are the definitions for some constants */ /******************************************************************************/ #define DIM 3 /* dimensionality of space */ #define SPECCODE 1 /* internal species code */ #define SPEC_NAME_LEN 64 /* max length of species name string */ #define SPEC_NAME_FMT "%63s" /* scan format of species name string */ /* Define prototype for Model Driver init */ int model_driver_create(KIM_ModelDriverCreate * const modelDriverCreate, KIM_LengthUnit const requestedLengthUnit, KIM_EnergyUnit const requestedEnergyUnit, KIM_ChargeUnit const requestedChargeUnit, KIM_TemperatureUnit const requestedTemperatureUnit, KIM_TimeUnit const requestedTimeUnit); /* Define prototypes for destroy */ /* defined as static to avoid namespace clashes with other codes */ static int destroy_routine(KIM_ModelDestroy * const modelDestroy); /* Define prototype for routines */ static int compute_routine(KIM_ModelCompute const * const modelCompute, KIM_ModelComputeArguments const * const modelComputeArguments); static int compute_arguments_create( KIM_ModelCompute const * const modelCompute, KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate); static int compute_arguments_destroy( KIM_ModelCompute const * const modelCompute, KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy); static int refresh_routine(KIM_ModelRefresh * const modelRefresh); static int write_parameterized_model(KIM_ModelWriteParameterizedModel const * const modelWriteParameterizedModel); /* Define model_buffer structure */ struct model_buffer { double influenceDistance; double cutoff; double cut_sq; int modelWillNotRequestNeighborsOfNoncontributingParticles; char species_name[SPEC_NAME_LEN]; double * params; }; /* compute function */ #undef KIM_LOGGER_FUNCTION_NAME #define KIM_LOGGER_FUNCTION_NAME KIM_ModelCompute_LogEntry #undef KIM_LOGGER_OBJECT_NAME #define KIM_LOGGER_OBJECT_NAME modelCompute static int compute_routine(KIM_ModelCompute const * const modelCompute, KIM_ModelComputeArguments const * const modelComputeArguments) { /* local variables */ double Rij; double Rik; double Rjk; double Rsqij; double Rsqik; double Rsqjk; double phi2; double dphi2_dRij; double phi2_contrib_prefactor; double phi2_force_contrib; double phi3; double dphi3_dRij; double dphi3_dRik; double dphi3_dRjk; double phi3_force_contrib_ij; double phi3_force_contrib_ik; double phi3_force_contrib_jk; double r_ij[DIM]; double r_ik[DIM]; double r_jk[DIM]; int ier; int atom_i; int neigh_count; int neigh_count2; int atom_j; int atom_k; int d; int const * neigh_list_of_current_part; struct model_buffer const * buffer; int comp_energy; int comp_force; int const * n_parts; int const * particle_species_codes; int const * particle_contributing; double const * cut_sq; double const * params; double const * coords; double * energy; double * force; int numb_of_part_neigh; /* get buffer from KIM object */ KIM_ModelCompute_GetModelBufferPointer(modelCompute, (void **) &buffer); /* unpack info from the buffer */ cut_sq = &(buffer->cut_sq); params = buffer->params; ier = KIM_ModelComputeArguments_GetArgumentPointerInteger( modelComputeArguments, KIM_COMPUTE_ARGUMENT_NAME_numberOfParticles, (int **) &n_parts) || KIM_ModelComputeArguments_GetArgumentPointerInteger( modelComputeArguments, KIM_COMPUTE_ARGUMENT_NAME_particleSpeciesCodes, (int **) &particle_species_codes) || KIM_ModelComputeArguments_GetArgumentPointerInteger( modelComputeArguments, KIM_COMPUTE_ARGUMENT_NAME_particleContributing, (int **) &particle_contributing) || KIM_ModelComputeArguments_GetArgumentPointerDouble( modelComputeArguments, KIM_COMPUTE_ARGUMENT_NAME_coordinates, (double **) &coords) || KIM_ModelComputeArguments_GetArgumentPointerDouble( modelComputeArguments, KIM_COMPUTE_ARGUMENT_NAME_partialEnergy, (double **) &energy) || KIM_ModelComputeArguments_GetArgumentPointerDouble( modelComputeArguments, KIM_COMPUTE_ARGUMENT_NAME_partialForces, (double **) &force); if (ier) { LOG_ERROR("Unable to get argument pointer."); return ier; } comp_energy = (energy != NULL); comp_force = (force != NULL); /* Check to be sure that the species are correct */ /**/ ier = TRUE; /* assume an error */ for (atom_i = 0; atom_i < *n_parts; ++atom_i) { if (SPECCODE != particle_species_codes[atom_i]) { LOG_ERROR("Unexpected species code detected."); return ier; } } ier = FALSE; /* everything is ok */ /* initialize potential energies, forces, and virial term */ if (comp_energy) { *energy = 0.0; } if (comp_force) { for (atom_i = 0; atom_i < *n_parts; ++atom_i) { for (d = 0; d < DIM; ++d) { force[atom_i * DIM + d] = 0.0; } } } /* Compute energy and forces */ /* loop over particles and compute energy and forces */ for (atom_i = 0; atom_i < *n_parts; ++atom_i) { if (particle_contributing[atom_i]) { ier = KIM_ModelComputeArguments_GetNeighborList( modelComputeArguments, 0, atom_i, &numb_of_part_neigh, &neigh_list_of_current_part); if (ier) { /* some sort of problem, exit */ LOG_ERROR("Unable to get neighbor list."); ier = TRUE; return ier; } /* loop over the neighbors of particle atom_i */ for (neigh_count = 0; neigh_count < numb_of_part_neigh; ++neigh_count) { atom_j = neigh_list_of_current_part[neigh_count]; /* get neighbor ID */ if (atom_j == atom_i) continue; /* compute relative position vector between atoms i and j, and squared * distance */ Rsqij = 0.0; for (d = 0; d < DIM; ++d) { r_ij[d] = coords[atom_j * DIM + d] - coords[atom_i * DIM + d]; /* compute squared distance */ Rsqij += r_ij[d] * r_ij[d]; } /* compute energy and force */ if (Rsqij < *cut_sq) { /* particles i and j are interacting */ Rij = sqrt(Rsqij); if (!(particle_contributing[atom_j] && atom_j < atom_i)) /* effective half-list */ { /* Account for whether this neighbor is a ghost atom (which would not * contribute any energy) */ phi2_contrib_prefactor = (particle_contributing[atom_j]) ? 1.0 : 0.5; if (comp_force) { /* compute two-body energy and its derivative */ calc_phi2_dphi2(params, Rij, &phi2, &dphi2_dRij); } else { /* compute just two-body energy */ calc_phi2_dphi2(params, Rij, &phi2, NULL); } /* contribution to energy */ if (comp_energy) { *energy += phi2_contrib_prefactor * phi2; } /* contribution to forces */ if (comp_force) { for (d = 0; d < DIM; ++d) { phi2_force_contrib = phi2_contrib_prefactor * dphi2_dRij * r_ij[d] / Rij; force[atom_i * DIM + d] += phi2_force_contrib; /* accumulate force on atom_i */ force[atom_j * DIM + d] -= phi2_force_contrib; /* accumulate force on atom_j */ } } } /* end effective half-list check (i < j) */ /* Secondary loop over neighbors to calculate three-body * interactions */ for (neigh_count2 = neigh_count + 1; neigh_count2 < numb_of_part_neigh; ++neigh_count2) { atom_k = neigh_list_of_current_part[neigh_count2]; /* get neighbor ID */ if (atom_k == atom_i) continue; /* compute relative position vector between atoms i and k, and * squared distance */ Rsqik = 0.0; for (d = 0; d < DIM; ++d) { r_ik[d] = coords[atom_k * DIM + d] - coords[atom_i * DIM + d]; /* compute squared distance */ Rsqik += r_ik[d] * r_ik[d]; } /* compute energy and force */ if (Rsqik < *cut_sq) { /* particles i and k are interacting */ Rik = sqrt(Rsqik); /* compute relative position vector between atoms j and k, and * squared distance */ Rsqjk = 0.0; for (d = 0; d < DIM; ++d) { r_jk[d] = coords[atom_k * DIM + d] - coords[atom_j * DIM + d]; /* compute squared distance */ Rsqjk += r_jk[d] * r_jk[d]; } if (SKIP_CHECK_ON_RJK || Rsqjk < *cut_sq) { /* particles j and k are interacting */ Rjk = sqrt(Rsqjk); if (comp_force) { /* compute three-body energy and its derivatives with * respect to Rij, Rik, and Rjk */ calc_phi3_dphi3(params, Rij, Rik, Rjk, &phi3, &dphi3_dRij, &dphi3_dRik, &dphi3_dRjk); } else { /* compute just three-body energy */ calc_phi3_dphi3( params, Rij, Rik, Rjk, &phi3, NULL, NULL, NULL); } /* contribution to energy */ if (comp_energy) { *energy += phi3; } /* contribution to forces */ if (comp_force) { for (d = 0; d < DIM; ++d) { phi3_force_contrib_ij = dphi3_dRij * r_ij[d] / Rij; phi3_force_contrib_ik = dphi3_dRik * r_ik[d] / Rik; phi3_force_contrib_jk = dphi3_dRjk * r_jk[d] / Rjk; force[atom_i * DIM + d] += phi3_force_contrib_ij + phi3_force_contrib_ik; /* accumulate force on atom_i */ force[atom_j * DIM + d] += -phi3_force_contrib_ij + phi3_force_contrib_jk; /* accumulate force on atom_j */ force[atom_k * DIM + d] += -phi3_force_contrib_ik - phi3_force_contrib_jk; /* accumulate force on atom_k */ } } } /* Optional check on if Rsqjk < *cut_sq */ } /* if Rsqik < *cut_sq */ } /* loop on neigh_count2 */ } /* if Rsqij < *cut_sq */ } /* loop on neigh_count */ } /* if particle_contributing */ } /* infinite while loop (terminated by break statements above) */ /* everything is great */ ier = FALSE; return ier; } /* Create function */ #undef KIM_LOGGER_FUNCTION_NAME #define KIM_LOGGER_FUNCTION_NAME KIM_ModelDriverCreate_LogEntry #undef KIM_LOGGER_OBJECT_NAME #define KIM_LOGGER_OBJECT_NAME modelDriverCreate int model_driver_create(KIM_ModelDriverCreate * const modelDriverCreate, KIM_LengthUnit const requestedLengthUnit, KIM_EnergyUnit const requestedEnergyUnit, KIM_ChargeUnit const requestedChargeUnit, KIM_TemperatureUnit const requestedTemperatureUnit, KIM_TimeUnit const requestedTimeUnit) { /* KIM variables */ int number_of_parameter_files; char const * param_file_1_name; /* Local variables */ FILE * fid; char species_name_string[SPEC_NAME_LEN]; KIM_SpeciesName species_name; int ier; int param_count; double * params; double conversion_factor; struct model_buffer * buffer; /* Use function pointer definitions to verify prototypes */ KIM_ModelDriverCreateFunction * create = model_driver_create; KIM_ModelComputeArgumentsCreateFunction * ca_create = compute_arguments_create; KIM_ModelComputeFunction * compute = compute_routine; KIM_ModelRefreshFunction * refresh = refresh_routine; KIM_ModelWriteParameterizedModelFunction * writeModel = write_parameterized_model; KIM_ModelComputeArgumentsDestroyFunction * ca_destroy = compute_arguments_destroy; KIM_ModelDestroyFunction * destroy = destroy_routine; (void) create; /* avoid unused parameter warnings */ (void) requestedChargeUnit; (void) requestedTemperatureUnit; (void) requestedTimeUnit; /* using requested units */ ier = KIM_ModelDriverCreate_SetUnits(modelDriverCreate, requestedLengthUnit, requestedEnergyUnit, KIM_CHARGE_UNIT_unused, KIM_TEMPERATURE_UNIT_unused, KIM_TIME_UNIT_unused); if (ier == TRUE) { LOG_ERROR("Unable to set units."); return ier; } ier = KIM_ModelDriverCreate_SetModelNumbering(modelDriverCreate, KIM_NUMBERING_zeroBased); if (ier == TRUE) { LOG_ERROR("Unable to set numbering."); return ier; } /* store pointer to functions in KIM object */ ier = KIM_ModelDriverCreate_SetRoutinePointer( modelDriverCreate, KIM_MODEL_ROUTINE_NAME_ComputeArgumentsCreate, KIM_LANGUAGE_NAME_c, TRUE, (KIM_Function *) ca_create) || KIM_ModelDriverCreate_SetRoutinePointer( modelDriverCreate, KIM_MODEL_ROUTINE_NAME_Compute, KIM_LANGUAGE_NAME_c, TRUE, (KIM_Function *) compute) || KIM_ModelDriverCreate_SetRoutinePointer( modelDriverCreate, KIM_MODEL_ROUTINE_NAME_Refresh, KIM_LANGUAGE_NAME_c, TRUE, (KIM_Function *) refresh) || KIM_ModelDriverCreate_SetRoutinePointer( modelDriverCreate, KIM_MODEL_ROUTINE_NAME_WriteParameterizedModel, KIM_LANGUAGE_NAME_c, FALSE, (KIM_Function *) writeModel) || KIM_ModelDriverCreate_SetRoutinePointer( modelDriverCreate, KIM_MODEL_ROUTINE_NAME_ComputeArgumentsDestroy, KIM_LANGUAGE_NAME_c, TRUE, (KIM_Function *) ca_destroy) || KIM_ModelDriverCreate_SetRoutinePointer( modelDriverCreate, KIM_MODEL_ROUTINE_NAME_Destroy, KIM_LANGUAGE_NAME_c, TRUE, (KIM_Function *) destroy); if (ier == TRUE) { LOG_ERROR("Unable to register model function pointers."); return ier; } /* get number of parameter files */ KIM_ModelDriverCreate_GetNumberOfParameterFiles(modelDriverCreate, &number_of_parameter_files); /* set param_file_1_name */ if (number_of_parameter_files != 1) { ier = TRUE; LOG_ERROR("Incorrect number of parameter files."); return ier; } ier = KIM_ModelDriverCreate_GetParameterFileName( modelDriverCreate, 0, ¶m_file_1_name); if (ier == TRUE) { LOG_ERROR("Unable to get parameter file name."); return ier; } /* Read in model parameters from parameter file */ fid = fopen(param_file_1_name, "r"); if (fid == NULL) { ier = TRUE; LOG_ERROR("Unable to open parameter file."); return ier; } ier = fscanf(fid, SPEC_NAME_FMT, species_name_string); /* species symbol */ /* Check that we successfully read in this line */ if (1 != ier) { ier = TRUE; fclose(fid); LOG_ERROR("Unable to read species from parameter file."); } params = (double *) malloc(sizeof(double) * NUM_PARAMS); if (params == NULL) { ier = TRUE; fclose(fid); LOG_ERROR("Unable to malloc memory for parameters."); return ier; } for (param_count = 0; param_count < NUM_PARAMS; ++param_count) { ier = fscanf(fid, "%lf", ¶ms[param_count]); /* Check that we successfully read in this line */ if (1 != ier) { ier = TRUE; LOG_ERROR("Unable to read parameter file."); } /* convert units */ if ((param_units[param_count][0] != 0.0) || (param_units[param_count][1] != 0.0)) { ier = KIM_ModelDriverCreate_ConvertUnit(*length_unit, *energy_unit, KIM_CHARGE_UNIT_unused, KIM_TEMPERATURE_UNIT_unused, KIM_TIME_UNIT_unused, requestedLengthUnit, requestedEnergyUnit, KIM_CHARGE_UNIT_unused, KIM_TEMPERATURE_UNIT_unused, KIM_TIME_UNIT_unused, param_units[param_count][0], param_units[param_count][1], 0.0, 0.0, 0.0, &conversion_factor); if (ier == TRUE) { fclose(fid); LOG_ERROR("Unable to convert units of parameter."); return ier; } params[param_count] *= conversion_factor; } } fclose(fid); /* register species */ species_name = KIM_SpeciesName_FromString(species_name_string); ier = KIM_ModelDriverCreate_SetSpeciesCode( modelDriverCreate, species_name, SPECCODE); if (ier == TRUE) { LOG_ERROR("Unable to set species code."); return ier; } /* allocate buffer */ buffer = (struct model_buffer *) malloc(sizeof(struct model_buffer)); if (NULL == buffer) { ier = TRUE; LOG_ERROR("Unable to malloc memory for buffer."); return ier; } /* setup buffer */ buffer->params = params; buffer->influenceDistance = get_influence_distance(buffer->params); buffer->cutoff = buffer->influenceDistance; buffer->cut_sq = buffer->cutoff * buffer->cutoff; buffer->modelWillNotRequestNeighborsOfNoncontributingParticles = 1; sprintf(buffer->species_name, "%s", /* no buffer overrun possible, due to above code */ species_name_string); /* end setup buffer */ /* store in model buffer */ KIM_ModelDriverCreate_SetModelBufferPointer(modelDriverCreate, (void *) buffer); /* publish model parameters */ for (param_count = 0; param_count < NUM_PARAMS; ++param_count) { ier = KIM_ModelDriverCreate_SetParameterPointerDouble( modelDriverCreate, 1, &(buffer->params[param_count]), param_strings[param_count][0], param_strings[param_count][1]); if (ier == TRUE) { LOG_ERROR("Unable to set parameter pointer(s)."); return TRUE; } } /* store model cutoff in KIM object */ KIM_ModelDriverCreate_SetInfluenceDistancePointer( modelDriverCreate, &(buffer->influenceDistance)); KIM_ModelDriverCreate_SetNeighborListPointers( modelDriverCreate, 1, &(buffer->cutoff), &(buffer->modelWillNotRequestNeighborsOfNoncontributingParticles)); return FALSE; } /* Refresh function */ #undef KIM_LOGGER_FUNCTION_NAME #define KIM_LOGGER_FUNCTION_NAME KIM_ModelRefresh_LogEntry #undef KIM_LOGGER_OBJECT_NAME #define KIM_LOGGER_OBJECT_NAME modelRefresh int refresh_routine(KIM_ModelRefresh * const modelRefresh) { struct model_buffer * buffer; /* get model buffer from KIM object */ KIM_ModelRefresh_GetModelBufferPointer(modelRefresh, (void **) &buffer); /* update cutoff and cutoff square */ buffer->influenceDistance = get_influence_distance(buffer->params); buffer->cutoff = buffer->influenceDistance; buffer->cut_sq = buffer->cutoff * buffer->cutoff; /* store model cutoff in KIM object */ KIM_ModelRefresh_SetInfluenceDistancePointer(modelRefresh, &(buffer->influenceDistance)); KIM_ModelRefresh_SetNeighborListPointers( modelRefresh, 1, &(buffer->cutoff), &(buffer->modelWillNotRequestNeighborsOfNoncontributingParticles)); return FALSE; } /* destroy function */ static int destroy_routine(KIM_ModelDestroy * const modelDestroy) { /* Local variables */ struct model_buffer * buffer; int ier; /* get model buffer from KIM object */ KIM_ModelDestroy_GetModelBufferPointer(modelDestroy, (void **) &buffer); /* free the buffer */ free(buffer->params); free(buffer); ier = FALSE; return ier; } /* compute arguments create routine */ #undef KIM_LOGGER_FUNCTION_NAME #define KIM_LOGGER_FUNCTION_NAME KIM_ModelComputeArgumentsCreate_LogEntry #undef KIM_LOGGER_OBJECT_NAME #define KIM_LOGGER_OBJECT_NAME modelComputeArgumentsCreate static int compute_arguments_create( KIM_ModelCompute const * const modelCompute, KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate) { int ier; (void) modelCompute; /* avoid unused parameter warning */ /* register arguments */ ier = KIM_ModelComputeArgumentsCreate_SetArgumentSupportStatus( modelComputeArgumentsCreate, KIM_COMPUTE_ARGUMENT_NAME_partialEnergy, KIM_SUPPORT_STATUS_optional) || KIM_ModelComputeArgumentsCreate_SetArgumentSupportStatus( modelComputeArgumentsCreate, KIM_COMPUTE_ARGUMENT_NAME_partialParticleEnergy, KIM_SUPPORT_STATUS_optional) || KIM_ModelComputeArgumentsCreate_SetArgumentSupportStatus( modelComputeArgumentsCreate, KIM_COMPUTE_ARGUMENT_NAME_partialForces, KIM_SUPPORT_STATUS_optional); if (ier == TRUE) { LOG_ERROR("Unable to set argument supportStatus."); return TRUE; } else { return FALSE; } } /* compute arguments destroy routine */ static int compute_arguments_destroy( KIM_ModelCompute const * const modelCompute, KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy) { (void) modelCompute; /* avoid unused parameter warning */ (void) modelComputeArgumentsDestroy; /* Nothing further to do */ return FALSE; } /* write parameterized model routine */ #undef KIM_LOGGER_FUNCTION_NAME #define KIM_LOGGER_FUNCTION_NAME KIM_ModelWriteParameterizedModel_LogEntry #undef KIM_LOGGER_OBJECT_NAME #define KIM_LOGGER_OBJECT_NAME modelWriteParameterizedModel /* */ #define STR_LEN 2048 static int write_parameterized_model( KIM_ModelWriteParameterizedModel const * const modelWriteParameterizedModel) { FILE * fp; char string_buffer[STR_LEN]; struct model_buffer const * buffer; char const * path; char const * model_name; int param_count; int max_str_len; /* get buffer from KIM object */ KIM_ModelWriteParameterizedModel_GetModelBufferPointer( modelWriteParameterizedModel, (void **) &buffer); KIM_ModelWriteParameterizedModel_GetPath(modelWriteParameterizedModel, &path); KIM_ModelWriteParameterizedModel_GetModelName(modelWriteParameterizedModel, &model_name); max_str_len = strlen(path) + strlen(model_name) + 9; if (max_str_len >= STR_LEN) { LOG_ERROR("Path and model name too long for internal buffers."); return TRUE; } sprintf(string_buffer, "%s.params", model_name); KIM_ModelWriteParameterizedModel_SetParameterFileName( modelWriteParameterizedModel, string_buffer); sprintf(string_buffer, "%s/%s.params", path, model_name); fp = fopen(string_buffer, "w"); if (NULL == fp) { LOG_ERROR("Unable to open parameter file for write."); return TRUE; } fprintf(fp, "%s\n", buffer->species_name); for (param_count = 0; param_count < NUM_PARAMS; ++param_count) { fprintf(fp, "%20.15f\n", buffer->params[param_count]); } fclose(fp); return FALSE; }