arcanis/moldyn
1
0
Fork 0
mirror of https://github.com/arcan1s/moldyn.git synced 2025-07-16 07:10:00 +00:00
Code Issues Projects Releases Wiki Activity

Added prefix 'mm_'

Browse Source
This commit is contained in:
arcan1s
2013-08-03 04:41:20 +04:00
parent e9e0e82c4d
commit ae379fd22d
96 changed files with 44 additions and 44 deletions
Download Patch File Download Diff File Expand all files Collapse all files

29
mm_radf/src/CMakeLists.txt Normal file
View File

@ -0,0 +1,29 @@
set ("${PROJECT}_VERSION_MAJOR" 1)
set ("${PROJECT}_VERSION_MINOR" 0)
set ("${PROJECT}_VERSION_PATCH" 1)
set ("${PROJECT}_VERSION" ${${PROJECT}_VERSION_MAJOR}.${${PROJECT}_VERSION_MINOR}.${${PROJECT}_VERSION_PATCH})
message (STATUS "${PROJECT}: Version ${${PROJECT}_VERSION}")
# set files
aux_source_directory (. SOURCES)
# set library
if (CMAKE_COMPILER_IS_GNUCXX)
set (ADDITIONAL_LIB m)
else ()
set (ADDITIONAL_LIB)
endif()
# message
message (STATUS "SOURCES: ${SOURCES}")
# link libraries and compile
add_executable (${PROJECT} ${SOURCES})
target_link_libraries (${PROJECT} ${ADDITIONAL_LIB})
# install properties
INSTALL (TARGETS ${PROJECT} DESTINATION bin)
if (ADD_INCLUDE)
INSTALL (FILES ${PUBLIC_HEADERS} DESTINATION include/${PROJECT})
endif ()

21
mm_radf/src/Makefile Normal file
View File

@ -0,0 +1,21 @@
PROJECT=RADF
CC=gcc
CFLAGS=-c -Wall -fPIC
LDFLAGS=-lm
SOURCES_DIR=src
SOURCES=main.c add_main.c coords.c messages.c radf.c radf_proc.c
OBJECTS=$(SOURCES:.c=.o)
EXECUTABLE=radf
$(PROJECT): $(SOURCES) $(EXECUTABLE)
$(EXECUTABLE): $(OBJECTS)
$(CC) $(LDFLAGS) $(OBJECTS) -o $@
.c.o:
$(CC) $(CFLAGS) $< -o $@
clean:
rm -f *.o

211
mm_radf/src/add_main.c Normal file
View File

@ -0,0 +1,211 @@
/**
* @file
*/
#include <stdio.h>
#include "messages.h"
/**
* @fn error_checking
*/
int error_checking (const float *cell, const int from, const char *input,
const int num_needed_at, const int *needed_at,
const char *output, const int to)
/**
* @brief function that checks errors in input variables
* @code
* error_checking (cell, from, input, num_needed_at, needed_at, output, to);
* @endcode
*
* @param cell massive of cell size
* @param from first trajectory step
* @param input input file name
* @param num_needed_at number of needed atom types
* @param needed_at massive of number of needed atom types
* @param output output file name
* @param to last trajectory step
*
* @return 11 - error in 'cell'
* @return 12 - error in 'input'
* @return 13 - error in 'output'
* @return 14 - error in 'from' or 'to'
* @return 15 - error in 'num_needed_at'
* @return 16 - error in 'needed_at'
* @return 0 - exit without errors
*/
{
if ((cell[0] == 0.0) || (cell[1] == 0.0) || (cell[2] == 0.0))
return 11;
if (input[0] == '#')
return 12;
if (output[0] == '#')
return 13;
if ((from == -1) || (to == -1))
return 14;
if ((num_needed_at != 2) && (num_needed_at != 6))
return 15;
if (num_needed_at == 6)
if ((needed_at[0] == needed_at[1]) ||
(needed_at[0] == needed_at[2]) ||
(needed_at[1] == needed_at[2]) ||
(needed_at[3] == needed_at[4]) ||
(needed_at[3] == needed_at[5]) ||
(needed_at[4] == needed_at[5]))
return 16;
return 0;
}
/**
* @fn print_message
*/
int print_message (const int quiet, FILE *std_output, const int log, FILE *f_log,
const int mode, const char *str)
/**
* @brief function that prints message in log and stdout
* @code
* print_message (quiet, stdout, log, f_log, 0, str);
* @endcode
*
* @param quiet status of quiet-mode
* @param std_output stdout
* @param log status of log-mode
* @param f_log log file
* @param mode number of message in "messages.c"
* @param str additional text in message
*
* @return 0 - exit without errors
*/
{
if ((quiet != 1) && (std_output != stderr))
message (0, mode, str, std_output);
if (log == 1)
message (1, mode, str, f_log);
return 0;
}
/**
* @fn printing_head
*/
int printing_head (const char *output, const int log, const int quiet, const int matrix,
const char *input, const int from, const int to, const float *cell,
const int mode, const double r_min, const double r_max,
const double r_step, const double ang_min, const double ang_max,
const double ang_step, const int *needed_at)
/**
* @brief function that prints header in output file
* @code
* printing_head (output, log, quiet, matrix, input, from, to, cell, mode, r_min,
* r_max, r_step, ang_min, ang_max, ang_step, needed_at);
* @endcode
*
* @param output output file nams
* @param log status of log-mode
* @param quiet status of quiet-mode
* @param matrix status of matrix-mode
* @param input mask of trajectory files
* @param from first trajectory step
* @param to last trajectory step
* @param cell massive of cell size
* @param mode 0 - if RDF, 1 - if RDF for center mass, 2 - if RADF
* @param r_max maximal radius
* @param r_min minimal radius
* @param r_step radius step
* @param ang_max maximal angle for RADF
* @param ang_min minimal angle for RADF
* @param ang_step anlge step for RADF
* @param needed_at massive of number of needed atom types
*
* @return 0 - exit without errors
*/
{
FILE *f_out;
f_out = fopen (output, "w");
fprintf (f_out, "radf ::: V.1.0.1 ::: 2013-07-23\n\n");
fprintf (f_out, "CONFIGURATION\n");
fprintf (f_out, "LOG=%i\nQUIET=%i\nMATRIX=%i\n", log, quiet, matrix);
fprintf (f_out, "MASK=%s\nFIRST=%i\nLAST=%i\n", input, from, to);
fprintf (f_out, "CELL=%.4f,%.4f,%.4f\n", cell[0], cell[1], cell[2]);
fprintf (f_out, "MODE=%i\n", mode);
fprintf (f_out, "R_MIN=%.3f\nR_MAX=%.3f\nR_STEP=%.3f", r_min, r_max, r_step);
switch (mode)
{
case 0:
fprintf (f_out, "ATOM=%i-%i\n", needed_at[0], needed_at[1]);
break;
case 1:
fprintf (f_out, "ATOM=%i,%i,%i-%i,%i,%i\n", needed_at[0], needed_at[1],
needed_at[2], needed_at[3], needed_at[4], needed_at[5]);
break;
case 2:
fprintf (f_out, "ANG_MIN=%.2f\nANG_MAX=%.2f\nANG_STEP=%.2f\n\
ATOM=%i,%i,%i-%i,%i,%i\n",ang_min, ang_max, ang_step, needed_at[0], needed_at[1],
needed_at[2], needed_at[3], needed_at[4], needed_at[5]);
break;
}
fprintf (f_out, "END\n\n");
fclose (f_out);
return 0;
}
/**
* @fn set_defaults
*/
int set_defaults (float *ang_max, float *ang_min, float *ang_step, float *cell,
int *from, char *input, int *log, int *matrix, float *r_max,
float *r_min, float *r_step, char *output, int *to, int *quiet)
/**
* @brief function that sets default values of variables
* @code
* set_defaults (&ang_max, &ang_min, &ang_step, cell, &from, input, &log, &r_max,
* &r_min, &r_step, output, &to, &quiet);
* @endcode
*
* @param ang_max maximal angle for RADF
* @param ang_min minimal angle for RADF
* @param ang_step angle step
* @param cell massive of cell size
* @param from first trajectory step
* @param input mask of trajectory files
* @param log status of log-mode
* @param matrix status of matrix-mode
* @param r_max maximal radius
* @param r_min minimal radius
* @param r_step radius step
* @param output output file name
* @param to last trajectory step
* @param quiet status of quiet-mode
*
* @return 0 - exit without errors
*/
{
int i;
*ang_max = 90.0;
*ang_min = 0.0;
*ang_step = 0.0;
for (i=0; i<3; i++)
cell[i] = 0.0;
*from = -1;
input[0] = '#';
*log = 0;
*matrix = 0;
*r_max = 15.0;
*r_min = 2.0;
*r_step = 0.2;
output[0] = '#';
*to = -1;
*quiet = 0;
return 0;
}

31
mm_radf/src/add_main.h Normal file
View File

@ -0,0 +1,31 @@
/**
* @file
*/
#ifndef ADD_MAIN_H
#define ADD_MAIN_H
/**
* @fn error_checking
*/
/**
* @fn print_message
*/
/**
* @fn printing_head
*/
/**
* @fn set_defaults
*/
int error_checking (const float *, const int, const char *, const int, const int *,
const char *, const int);
int print_message (const int, FILE *, const int, FILE *, const int, const char *);
int printing_head (const char *, const int, const int, const int, const char *,
const int, const int, const float *, const int, const double,
const double, const double, const double, const double, const double,
const int *);
int set_defaults (float *, float *, float *, float *, int *, char *, int *, int *,
float *, float *, float *, char *, int *, int *);
#endif /* ADD_MAIN_H */

354
mm_radf/src/coords.c Normal file
View File

@ -0,0 +1,354 @@
/**
* @file
*/
#include <stdio.h>
#include <stdlib.h>
/**
* @fn reading_coords
*/
int reading_coords (const int mode, const char *filename, const int type_inter,
const int *label_atom, const float *cell, int *num_mol,
int *num_atoms, int *true_label_mol, int *label_mol,
int *type_atoms, float *coords, char *ch_type_atoms)
/**
* @brief function that reads coordinates from special file format
* @code
* reading_coords (0, filename, type_inter, label_atom, cell, &num_mol, &num_atoms,
* true_label_mol, label_mol, type_atoms, coords, ch_type_atoms);
* @endcode
*
* @param mode mode of reading; '1' is statgen, '2' is envir or
* frad, '3' is agl
* @param filename input file name
* @param type_inter number of needed atoms
* (number of needed molecules)
* @param label_atom massive of needed atom types
* (massive of needed molecules)
* @param cell massive of cell size
* @param num_mol number of molecules
* @param num_atoms number of atoms
* @param true_label_mol massive of true numbers of molecule for atoms
* @param label_mol massive of numbers of molecule for atoms
* @param type_atoms massive of atom types
* @param coords massive of coordinates
* @param ch_type_atoms massive of char atom types
*
* @return 1 - file $filename does not exist
* @return 2 - unknown mode
* @return 0 - exit without errors
*/
{
char at_symb[32], file_string[256];
int atoms, cur_at_num, cur_at_type, cur_mol, i, j, tr_num_atoms, ref_mol, x, y;
float cur_coords[3], *not_tr_coords, ref[3];
FILE *inp;
/* cur_* temp variables
* at_symb temp variable
* file_string temp string variable
* atoms total number of atoms in system
* tr_num_atoms number of translated atoms (must be 8*num_atoms)
* ref_mol number of molecule for reference in translation
* not_tr_coords massive of not translated coordinates
* ref massive of coordinates of reference molecule
* inp input file
*/
/// <b>Work blocks</b>
*num_atoms = 0;
*num_mol = 0;
/// <pre> reading file </pre>
inp = fopen (filename, "r");
if (inp == NULL)
return 1;
ref_mol = -1;
fscanf (inp, "%i", &atoms);
not_tr_coords = (float *) malloc (3 * atoms * sizeof (float));
fgets (file_string, 256, inp);
for (i=0; i<atoms; i++)
{
fgets (file_string, 256, inp);
sscanf (file_string, "%i%s%f%f%f%i%i", &cur_at_num, at_symb, &cur_coords[0],
&cur_coords[1], &cur_coords[2], &cur_at_type, &cur_mol);
// reading variables according to selected mode
switch (mode)
{
case 0:
// mode == 0 (selected atoms)
for (j=0; j<type_inter; j++)
if (cur_at_type == label_atom[j])
{
not_tr_coords[3**num_atoms+0] = cur_coords[0];
not_tr_coords[3**num_atoms+1] = cur_coords[1];
not_tr_coords[3**num_atoms+2] = cur_coords[2];
if (ref_mol != cur_mol)
{
ref_mol = cur_mol;
true_label_mol[*num_mol] = ref_mol;
*num_mol = *num_mol + 1;
}
label_mol[*num_atoms] = *num_mol - 1;
type_atoms[*num_atoms] = j;
*num_atoms = *num_atoms + 1;
}
break;
case 1:
// mode == 1 (all atoms)
not_tr_coords[3**num_atoms+0] = cur_coords[0];
not_tr_coords[3**num_atoms+1] = cur_coords[1];
not_tr_coords[3**num_atoms+2] = cur_coords[2];
ch_type_atoms[2**num_atoms+0] = at_symb[0];
ch_type_atoms[2**num_atoms+1] = at_symb[1];
if (ref_mol != cur_mol)
{
ref_mol = cur_mol;
true_label_mol[*num_mol] = ref_mol;
*num_mol = *num_mol + 1;
}
label_mol[*num_atoms] = *num_mol - 1;
type_atoms[*num_atoms] = j;
*num_atoms = *num_atoms + 1;
break;
case 2:
// mode == 2 (selected molecules)
for (j=0; j<type_inter; j++)
if (cur_mol == label_atom[j])
{
not_tr_coords[3**num_atoms+0] = cur_coords[0];
not_tr_coords[3**num_atoms+1] = cur_coords[1];
not_tr_coords[3**num_atoms+2] = cur_coords[2];
ch_type_atoms[2**num_atoms+0] = at_symb[0];
ch_type_atoms[2**num_atoms+1] = at_symb[1];
if (ref_mol != cur_mol)
{
ref_mol = cur_mol;
true_label_mol[*num_mol] = ref_mol;
*num_mol = *num_mol + 1;
}
label_mol[*num_atoms] = *num_mol - 1;
type_atoms[*num_atoms] = j;
*num_atoms = *num_atoms + 1;
}
break;
default: return 2;
}
}
fclose (inp);
/// <pre> translation </pre>
tr_num_atoms = *num_atoms;
for (i=0; i<*num_atoms; i++)
for (j=0; j<3; j++)
coords[3*i+j] = not_tr_coords[3*i+j];
// assign initial value to reference coordinates
ref_mol = label_mol[0];
for (i=0; i<3; i++)
ref[i] = coords[3*0+i];
for (i=0; i<*num_atoms; i++)
{
if (label_mol[i] != ref_mol)
{
ref_mol = label_mol[i];
for (j=0; j<3; j++)
ref[j] = not_tr_coords[3*i+j];
}
for (x=0; x<3; x++)
{
if (ref[x] >= 0.0)
// if xyz >= 0.0 A
{
for (j=0; j<3; j++)
if (j == x)
coords[3*tr_num_atoms+j] = not_tr_coords[3*i+j] - cell[j];
else
coords[3*tr_num_atoms+j] = not_tr_coords[3*i+j];
label_mol[tr_num_atoms] = label_mol[i];
type_atoms[tr_num_atoms] = type_atoms[i];
tr_num_atoms++;
}
else
// if xyz < 0.0 A
{
for (j=0; j<3; j++)
if (j == x)
coords[3*tr_num_atoms+j] = not_tr_coords[3*i+j] + cell[j];
else
coords[3*tr_num_atoms+j] = not_tr_coords[3*i+j];
label_mol[tr_num_atoms] = label_mol[i];
type_atoms[tr_num_atoms] = type_atoms[i];
tr_num_atoms++;
}
}
for (x=0; x<3; x++)
{
for (y=x+1; y<3; y++)
{
if ((ref[x] >= 0.0) && (ref[y] >= 0.0))
// if xyz and xyz >= 0.0 A
{
for (j=0; j<3; j++)
if ((j == x) || (j == y))
coords[3*tr_num_atoms+j] = not_tr_coords[3*i+j] - cell[j];
else
coords[3*tr_num_atoms+j] = not_tr_coords[3*i+j];
label_mol[tr_num_atoms] = label_mol[i];
type_atoms[tr_num_atoms] = type_atoms[i];
tr_num_atoms++;
}
if ((ref[x] < 0.0) && (ref[y] < 0.0))
// if xyz and xyz < 0.0 A
{
for (j=0; j<3; j++)
if ((j == x) || (j == y))
coords[3*tr_num_atoms+j] = not_tr_coords[3*i+j] + cell[j];
else
coords[3*tr_num_atoms+j] = not_tr_coords[3*i+j];
label_mol[tr_num_atoms] = label_mol[i];
type_atoms[tr_num_atoms] = type_atoms[i];
tr_num_atoms++;
}
}
for (y=0; y<3; y++)
if ((ref[x] < 0.0) && (ref[y] >= 0.0))
// if xyz OR xyz >= 0.0
{
for (j=0; j<3; j++)
{
if (j == x)
coords[3*tr_num_atoms+j] = not_tr_coords[3*i+j] + cell[j];
if (j == y)
coords[3*tr_num_atoms+j] = not_tr_coords[3*i+j] - cell[j];
if ((j != x) && (j != y))
coords[3*tr_num_atoms+j] = not_tr_coords[3*i+j];
}
label_mol[tr_num_atoms] = label_mol[i];
type_atoms[tr_num_atoms] = type_atoms[i];
tr_num_atoms++;
}
}
if ((ref[0] >= 0.0) && (ref[1] >= 0.0) && (ref[2] >= 0.0))
// if x and y and z >= 0.0 A
{
coords[3*tr_num_atoms+0] = not_tr_coords[3*i+0] - cell[0];
coords[3*tr_num_atoms+1] = not_tr_coords[3*i+1] - cell[1];
coords[3*tr_num_atoms+2] = not_tr_coords[3*i+2] - cell[2];
label_mol[tr_num_atoms] = label_mol[i];
type_atoms[tr_num_atoms] = type_atoms[i];
tr_num_atoms++;
}
if ((ref[0] >= 0.0) && (ref[1] >= 0.0) && (ref[2] < 0.0))
// if x and y >= 0.0 A and z < 0.0 A
{
coords[3*tr_num_atoms+0] = not_tr_coords[3*i+0] - cell[0];
coords[3*tr_num_atoms+1] = not_tr_coords[3*i+1] - cell[1];
coords[3*tr_num_atoms+2] = not_tr_coords[3*i+2] + cell[2];
label_mol[tr_num_atoms] = label_mol[i];
type_atoms[tr_num_atoms] = type_atoms[i];
tr_num_atoms++;
}
if ((ref[0] >= 0.0) && (ref[1] < 0.0) && (ref[2] >= 0.0))
// if x and z >= 0.0 A and y < 0.0 A
{
coords[3*tr_num_atoms+0] = not_tr_coords[3*i+0] - cell[0];
coords[3*tr_num_atoms+1] = not_tr_coords[3*i+1] + cell[1];
coords[3*tr_num_atoms+2] = not_tr_coords[3*i+2] - cell[2];
label_mol[tr_num_atoms] = label_mol[i];
type_atoms[tr_num_atoms] = type_atoms[i];
tr_num_atoms++;
}
if ((ref[0] < 0.0) && (ref[1] >= 0.0) && (ref[2] >= 0.0))
// if y and z >= 0.0 A and x < 0.0 A
{
coords[3*tr_num_atoms+0] = not_tr_coords[3*i+0] + cell[0];
coords[3*tr_num_atoms+1] = not_tr_coords[3*i+1] - cell[1];
coords[3*tr_num_atoms+2] = not_tr_coords[3*i+2] - cell[2];
label_mol[tr_num_atoms] = label_mol[i];
type_atoms[tr_num_atoms] = type_atoms[i];
tr_num_atoms++;
}
if ((ref[0] < 0.0) && (ref[1] < 0.0) && (ref[2] >= 0.0))
// if x and y < 0.0 A and z >= 0.0 A
{
coords[3*tr_num_atoms+0] = not_tr_coords[3*i+0] + cell[0];
coords[3*tr_num_atoms+1] = not_tr_coords[3*i+1] + cell[1];
coords[3*tr_num_atoms+2] = not_tr_coords[3*i+2] - cell[2];
label_mol[tr_num_atoms] = label_mol[i];
type_atoms[tr_num_atoms] = type_atoms[i];
tr_num_atoms++;
}
if ((ref[0] < 0.0) && (ref[1] >= 0.0) && (ref[2] < 0.0))
// if x and z < 0.0 A and y >= 0.0 A
{
coords[3*tr_num_atoms+0] = not_tr_coords[3*i+0] + cell[0];
coords[3*tr_num_atoms+1] = not_tr_coords[3*i+1] - cell[1];
coords[3*tr_num_atoms+2] = not_tr_coords[3*i+2] + cell[2];
label_mol[tr_num_atoms] = label_mol[i];
type_atoms[tr_num_atoms] = type_atoms[i];
tr_num_atoms++;
}
if ((ref[0] >= 0.0) && (ref[1] < 0.0) && (ref[2] < 0.0))
// if x >= 0.0 A and y and z < 0.0 A
{
coords[3*tr_num_atoms+0] = not_tr_coords[3*i+0] - cell[0];
coords[3*tr_num_atoms+1] = not_tr_coords[3*i+1] + cell[1];
coords[3*tr_num_atoms+2] = not_tr_coords[3*i+2] + cell[2];
label_mol[tr_num_atoms] = label_mol[i];
type_atoms[tr_num_atoms] = type_atoms[i];
tr_num_atoms++;
}
if ((ref[0] < 0.0) && (ref[1] < 0.0) && (ref[2] < 0.0))
// if x and y and z < 0.0 A
{
coords[3*tr_num_atoms+0] = not_tr_coords[3*i+0] + cell[0];
coords[3*tr_num_atoms+1] = not_tr_coords[3*i+1] + cell[1];
coords[3*tr_num_atoms+2] = not_tr_coords[3*i+2] + cell[2];
label_mol[tr_num_atoms] = label_mol[i];
type_atoms[tr_num_atoms] = type_atoms[i];
tr_num_atoms++;
}
}
/// <pre> free memory </pre>
free (not_tr_coords);
return 0;
}

16
mm_radf/src/coords.h Normal file
View File

@ -0,0 +1,16 @@
/**
* @file
*/
#ifndef COORDS_H
#define COORDS_H
/**
* @fn reading_coords
*/
int reading_coords (const int, const char *, const int, const int *,
const float *, int *, int *, int *, int *, int *, float *,
char *);
#endif /* COORDS_H */

420
mm_radf/src/main.c Normal file
View File

@ -0,0 +1,420 @@
/**
* @file
*/
/**
* @mainpage mm_radf
* @image latex ./logo.png
*
* @section intro_sec Introduction
*
* <b>About this program</b>:
* <ul>
* <li>Program that calculates radial distribution function (RDF) or radial-angles
* distribution function
* </ul>
*
* <b>Developer</b>:
* <ul>
* <li>Evgeniy Alekseev aka arcanis <pre><esalexeev (at) gmail (dot) com></pre>
*</ul>
* <b>License</b>:
* <ul>
* <li>GPL
* </ul>
*
* @section How-To-Use How to use
* Usage:
* <pre>
* mm_radf -i INPUT -s FIRST,LAST -c X,Y,Z -a ... -o OUTPUT [ -r MIN,MAX ] [ -rs R_STEP ]
* [ -a MIN,MAX ] [ -as ANG_STEP ] [ -m ] [ -l LOGFILE ] [ -q ] [ -h ]
* Parametrs:
* -i - mask of input files
* -s - trajectory steps (integer)
* -c - cell size (float), A
* -a - atom types (integer). Format: 'ATOM1-ATOM2' or 'A1,A2,A3-B1,B2,B3'
* (will enable RDF calculation for center mass automaticaly)
* -o - output file name
* -r - minimal and maximal radii for analyze (float), A. Default is '2.0,15.0'
* -rs - radius step for analyze (float), A. Default is '0.2'
* -a - minimal and maximal angles for analyze (float), deg. Default is '0.0,90.0'
* -as - angle step for analyze (float), deg. This option will enable RADF
* calculation automaticaly
* -m - matrix output enable
* -l - log enable
* -q - quiet enable
* -h - show this help and exit
* </pre>
*
* @page Install
*
* @section Requirements Requirements
* The application mm_radf requires the following external stuff:
* - cmake >= 2.8
* - gcc >= 4.8
*
* @section How-To How to install
*
* @subsection Linux Linux
* @code
* mkdir build && cd build
* cmake -DCMAKE_INSTALL_PREFIX=/usr -DCMAKE_BUILD_TYPE=Release ../
* make
* make install
* @endcode
*
* @subsection Windows Windows
* @code
* create project file using 'cmake'
* compile project
* @endcode
* You may also download compiled executable file for Win_x86.
*
* @page Changelog
* V.1.0.1 (2013-07-27)
* * initial release
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "add_main.h"
#include "coords.h"
#include "messages.h"
#include "radf.h"
#include "radf_proc.h"
/**
* @fn main
*/
int main(int argc, char *argv[])
/**
* @return 1 - error in error_checking
* @return 2 - input file does not exist
* @return 3 - memory error
* @return 4 - unknown flag
* @return 0 - exit without errors
*/
{
char filename[256], tmp_str[2048];
int error, i, j;
FILE *f_inp, *f_log;
char input[256], logfile[256], output[256];
float ang_max, ang_min, ang_step, cell[3], *coords, r_max, r_min, r_step;
int from, *label_mol, log, matrix, mode, num_atoms, num_mol, num_needed_at,
needed_at[6], quiet, *radf, step, to, *true_label_mol, *type_atoms;
/* input input file name
* logfile log file name
* output output file name
*
* ang_max maximal angle for RADF
* ang_min minimal angle for RADF
* ang_step anlge step for RADF
* cell massive of cell size
* coords massive of coordinates
* r_max maximal radius
* r_min minimal radius
* r_step radius step
*
* from first trajectory step
* label_mol massive of numbers of molecule for atoms
* log status of log-mode
* matrix status of matrix-mode
* mode 0 - if RDF, 1 - if RDF for center mass, 2 - if RADF
* num_atoms number of atoms
* num_mol number of molecules
* num_needed_at number of needed atom types
* needed_at massive of number of needed atom types
* quiet status of quiet-mode
* radf not normed RADF
* step $(to - from + 1)
* to last trajectory step
* true_label_mol massive of true numbers of molecule for atoms
* type_atoms massive of atom types
*/
set_defaults (&ang_max, &ang_min, &ang_step, cell, &from, input, &log, &matrix,
&r_max, &r_min, &r_step, output, &to, &quiet);
for (i=1; i<argc; i++)
{
if ((argv[i][0] == '-') && (argv[i][1] == 'h') && (argv[i][2] == '\0'))
{
sprintf (tmp_str, " mm_radf\n");
sprintf (tmp_str, "%sProgram that calculates radial distribution function (RDF) or radial-angles\n", tmp_str);
sprintf (tmp_str, "%sdistribution function\n", tmp_str);
sprintf (tmp_str, "%sVersion : 1.0.1 License : GPL\n", tmp_str);
sprintf (tmp_str, "%s Evgeniy Alekseev aka arcanis\n", tmp_str);
sprintf (tmp_str, "%s E-mail : esalexeev@gmail.com\n\n", tmp_str);
sprintf (tmp_str, "%sUsage:\n", tmp_str);
sprintf (tmp_str, "%smm_radf -i INPUT -s FIRST,LAST -c X,Y,Z -at ... -o OUTPUT [ -r MIN,MAX ] [ -rs R_STEP ]\n", tmp_str);
sprintf (tmp_str, "%s [ -a MIN,MAX ] [ -as ANG_STEP ] [ -m ] [ -l LOGFILE ] [ -q ] [ -h ]\n\n", tmp_str);
sprintf (tmp_str, "%sParametrs:\n", tmp_str);
sprintf (tmp_str, "%s -i - mask of input files\n", tmp_str);
sprintf (tmp_str, "%s -s - trajectory steps (integer)\n", tmp_str);
sprintf (tmp_str, "%s -c - cell size (float), A\n", tmp_str);
sprintf (tmp_str, "%s -at - atom types (integer). Format: 'ATOM1-ATOM2' or 'A1,A2,A3-B1,B2,B3'\n", tmp_str);
sprintf (tmp_str, "%s (will enable RDF calculation for center mass automaticaly)\n", tmp_str);
sprintf (tmp_str, "%s -o - output file name\n", tmp_str);
sprintf (tmp_str, "%s -r - minimal and maximal radii for analyze (float), A. Default is '2.0,15.0'\n", tmp_str);
sprintf (tmp_str, "%s -rs - radius step for analyze (float), A. Default is '0.2'\n", tmp_str);
sprintf (tmp_str, "%s -a - minimal and maximal angles for analyze (float), deg. Default is '0.0,90.0'\n", tmp_str);
sprintf (tmp_str, "%s -as - angle step for analyze (float), deg. This option will enable RADF\n", tmp_str);
sprintf (tmp_str, "%s calculation automaticaly\n", tmp_str);
sprintf (tmp_str, "%s -m - matrix output enable\n", tmp_str);
sprintf (tmp_str, "%s -l - log enable\n", tmp_str);
sprintf (tmp_str, "%s -q - quiet enable\n", tmp_str);
sprintf (tmp_str, "%s -h - show this help and exit\n", tmp_str);
fputs (tmp_str, stdout);
return 0;
}
else if ((argv[i][0] == '-') && (argv[i][1] == 'i') && (argv[i][2] == '\0'))
// input file
{
strcpy (input, argv[i+1]);
i++;
}
else if ((argv[i][0] == '-') && (argv[i][1] == 's') && (argv[i][2] == '\0'))
// steps
{
sscanf (argv[i+1], "%i,%i", &from, &to);
if (from > to)
{
to += from;
from = to - from;
to -= from;
}
step = to - from + 1;
i++;
}
else if ((argv[i][0] == '-') && (argv[i][1] == 'c') && (argv[i][2] == '\0'))
// cell size
{
sscanf (argv[i+1], "%f,%f,%f", &cell[0], &cell[1], &cell[2]);
i++;
}
else if ((argv[i][0] == '-') && (argv[i][1] == 'a') && (argv[i][2] == 't') && (argv[i][3] == '\0'))
// atom types
{
num_needed_at = 2;
for (j=0; j<strlen(argv[i+1]); j++)
if (argv[i+1][j] == ',')
num_needed_at++;
if (num_needed_at == 2)
sscanf (argv[i+1], "%i-%i", &needed_at[0], &needed_at[1]);
else if (num_needed_at == 6)
sscanf (argv[i+1], "%i,%i,%i-%i,%i,%i", &needed_at[0], &needed_at[1], &needed_at[2],
&needed_at[3], &needed_at[4], &needed_at[5]);
i++;
}
else if ((argv[i][0] == '-') && (argv[i][1] == 'o') && (argv[i][2] == '\0'))
// output file
{
strcpy (output, argv[i+1]);
i++;
}
else if ((argv[i][0] == '-') && (argv[i][1] == 'r') && (argv[i][2] == '\0'))
// radii
{
sscanf (argv[i+1], "%f,%f", &r_min, &r_max);
if (r_min > r_max)
{
r_min += r_max;
r_max = r_min - r_max;
r_min -= r_max;
}
i++;
}
else if ((argv[i][0] == '-') && (argv[i][1] == 'r') && (argv[i][2] == 's') && (argv[i][3] == '\0'))
// radius step
{
sscanf (argv[i+1], "%f", &r_step);
i++;
}
else if ((argv[i][0] == '-') && (argv[i][1] == 'a') && (argv[i][2] == '\0'))
// angles
{
sscanf (argv[i+1], "%f,%f", &ang_min, &ang_max);
if (ang_min > ang_max)
{
ang_min += ang_max;
ang_max = ang_min - ang_max;
ang_min -= ang_max;
}
i++;
}
else if ((argv[i][0] == '-') && (argv[i][1] == 'a') && (argv[i][2] == 's') && (argv[i][3] == '\0'))
// angle step
{
sscanf (argv[i+1], "%f", &ang_step);
i++;
}
else if ((argv[i][0] == '-') && (argv[i][1] == 'm') && (argv[i][2] == '\0'))
// matrix mode
{
matrix = 1;
}
else if ((argv[i][0] == '-') && (argv[i][1] == 'l') && (argv[i][2] == '\0'))
// log mode
{
log = 1;
strcpy (logfile, argv[i+1]);
i++;
}
else if ((argv[i][0] == '-') && (argv[i][1] == 'q') && (argv[i][2] == '\0'))
// quiet mode
{
quiet = 1;
}
else
// unknown flag
{
return 4;
}
}
if (log == 1)
f_log = fopen (logfile, "w");
print_message (quiet, stdout, log, f_log, 0, argv[0]);
print_message (quiet, stdout, log, f_log, 1, argv[0]);
// error check
error = error_checking (cell, from, input, num_needed_at, needed_at, output, to);
if (error != 0)
{
print_message (quiet, stderr, log, f_log, 17, argv[0]);
return 1;
}
print_message (quiet, stdout, log, f_log, 2, argv[0]);
// processing
// initial variables
// set mode
if ((num_needed_at == 6) && (ang_step != 0.0))
mode = 2;
else if (num_needed_at == 6)
mode = 1;
else if (num_needed_at == 2)
mode = 0;
sprintf (filename, "%s.%03i", input, from);
print_message (quiet, stdout, log, f_log, 3, filename);
f_inp = fopen (filename, "r");
if (f_inp == NULL)
{
print_message (quiet, stderr, log, f_log, 18, filename);
return 2;
}
fscanf (f_inp, "%i", &num_atoms);
fclose (f_inp);
coords = (float *) malloc (2 * 3 * 8 * num_atoms * sizeof (float));
label_mol = (int *) malloc (2 * 8 * num_atoms * sizeof (int));
if (mode == 2)
{
i = (r_max - r_min) / r_step;
j = (ang_max - ang_min) / ang_step;
i *= j;
}
else
i = (r_max - r_min) / r_step;
radf = (int *) malloc (i * sizeof (int));
for (j=0; j<i; j++)
radf[j] = 0;
true_label_mol = (int *) malloc (num_atoms * sizeof (int));
type_atoms = (int *) malloc (2 * 8 * num_atoms * sizeof (int));
// error checking
if ((coords == NULL) ||
(label_mol == NULL) ||
(radf == NULL) ||
(true_label_mol == NULL) ||
(type_atoms == NULL))
{
print_message (quiet, stderr, log, f_log, 19, argv[0]);
return 3;
}
sprintf (tmp_str, "%6cOutput file: '%s';\n%6cLog: %i;\n%6cQuiet: %i;\n%6cMatrix mode: %i;\n\
%6cMask: %s;\n%6cFirst step: %i;\n%6cLast step: %i;\n%6cCell size: %.4f, %.4f, %.4f;\n%6cMode: %i;\n\
%6cR_MIN: %6.3f; R_MAX: %6.3f; R_STEP: %6.3f;\n", ' ', output, ' ', log, ' ', quiet, ' ', matrix,
' ', input, ' ', from, ' ', to, ' ', cell[0], cell[1], cell[2], ' ', mode, ' ', r_min, r_max, r_step);
switch (mode)
{
case 0:
sprintf (tmp_str, "%s%6cATOM TYPES: %i-%i\n", tmp_str, ' ', needed_at[0], needed_at[1]);
break;
case 1:
sprintf (tmp_str, "%s%6cATOM TYPES: %i,%i,%i-%i,%i,%i\n", tmp_str, ' ', needed_at[0], needed_at[1],
needed_at[2], needed_at[3], needed_at[4], needed_at[5]);
break;
case 2:
sprintf (tmp_str, "%s%6cANG_MIN: %6.2f; ANG_MAX: %6.2f; ANG_STEP: %6.2f;\n\
%6cATOM TYPES: %i,%i,%i-%i,%i,%i\n", tmp_str, ' ', ang_min, ang_max, ang_step, ' ', needed_at[0],
needed_at[1], needed_at[2], needed_at[3], needed_at[4], needed_at[5]);
break;
}
print_message (quiet, stdout, log, f_log, 5, tmp_str);
// head
printing_head (output, log, quiet, matrix, input, from, to, cell, mode, r_min,
r_max, r_step, ang_min, ang_max, ang_step, needed_at);
print_message (quiet, stdout, log, f_log, 6, argv[0]);
// main cycle
for (i=from; i<to+1; i++)
{
sprintf (filename, "%s.%03i", input, i);
print_message (quiet, stdout, log, f_log, 7, filename);
error = reading_coords (0, filename, num_needed_at, needed_at, cell, &num_mol,
&num_atoms, true_label_mol, label_mol, type_atoms,
coords, tmp_str);
if (error == 0)
{
sprintf (tmp_str, "%6cNumber of molecules: %i; %6cNumber of atoms: %i\n",
' ', num_mol, ' ', num_atoms);
print_message (quiet, stdout, log, f_log, 8, tmp_str);
error = 1;
switch (mode)
{
case 0:
error = search_rdf (num_atoms, type_atoms, label_mol, coords, r_min,
r_max, r_step, radf);
break;
case 1:
error = search_rdf_centr (num_atoms, type_atoms, label_mol, coords, r_min,
r_max, r_step, radf);
if (error == 0)
case 2:
error = search_radf (num_atoms, type_atoms, label_mol, coords, r_min, r_max,
r_step, ang_min, ang_max, ang_step, radf);
break;
}
if (error == 0)
print_message (quiet, stdout, log, f_log, 22, argv[0]);
}
}
print_message (quiet, stdout, log, f_log, 13, argv[0]);
print_message (quiet, stdout, log, f_log, 14, output);
// tail
print_result (output, matrix, mode, step, num_atoms, r_min, r_max, r_step, ang_min,
ang_max, ang_step, cell, radf);
print_message (quiet, stdout, log, f_log, 15, argv[0]);
// free memory
free (coords);
free (label_mol);
free (radf);
free (type_atoms);
free (true_label_mol);
print_message (quiet, stdout, log, f_log, 16, argv[0]);
if (log == 1)
fclose (f_log);
return 0;
}

119
mm_radf/src/messages.c Normal file
View File

@ -0,0 +1,119 @@
/**
* @file
*/
#include <stdio.h>
#include <time.h>
/**
* @fn message
*/
int message (const int log, const int mode, const char *text, FILE *output)
/**
* @brief function that prints messages to output
* @code
* message (log, mode, text, output);
* @endcode
*
* @param log equal to 1 if print to logfile
* @param mode number of message
* @param text additional text
* @param output output file (may be stdout)
*
* @return 1 - unknown mode
* @return 0 - exit without errors
*/
{
char out[4096];
if (log == 1)
{
char time_str[256];
time_t t = time (NULL);
struct tm* aTm = localtime (&t);
sprintf (time_str, "[%04d-%02d-%02d %02d:%02d:%02d] [%2i]: ", aTm->tm_year+1900,
aTm->tm_mon+1, aTm->tm_mday, aTm->tm_hour, aTm->tm_min, aTm->tm_sec, mode);
fputs (time_str, output);
}
switch (mode)
{
case 0:
sprintf (out, "Start program: '%s'\n", text);
break;
case 1:
sprintf (out, "Checking errors\n");
break;
case 2:
sprintf (out, "Errors are not detected\n");
break;
case 3:
sprintf (out, "Initialization of variables from file '%s'\n", text);
break;
case 4:
sprintf (out, "%6cAglomerate was selected successfully\n", ' ');
break;
case 5:
sprintf (out, "Initial parametrs: \n%s", text);
break;
case 6:
sprintf (out, "Processing\n");
break;
case 7:
sprintf (out, "Open file: '%s'\n", text);
break;
case 8:
sprintf (out, "%s", text);
break;
case 9:
sprintf (out, "%6cSize of variables was changed successfully\n", ' ');
break;
case 10:
sprintf (out, "%6cConnectivity matrix was created successfully\n", ' ');
break;
case 11:
sprintf (out, "%6cConnectivity matrix was processed successfully\n", ' ');
break;
case 12:
sprintf (out, "%6cResult was printed to file '%s' successfully\n", ' ', text);
break;
case 13:
sprintf (out, "End of processing\n");
break;
case 14:
sprintf (out, "Print result to file '%s'\n", text);
break;
case 15:
sprintf (out, "Free memory\n");
break;
case 16:
sprintf (out, "Exiting without errors\n");
break;
case 17:
sprintf (out, "Something wrong!\nSee '%s -h' for more details\n", text);
break;
case 18:
sprintf (out, "File '%s' not found\nError\n", text);
break;
case 19:
sprintf (out, "Memory error\n");
break;
case 20:
sprintf (out, "%6cCenter of molecules was set successfully\n", ' ');
break;
case 21:
sprintf (out, "%6cEnvironment was selected successfully\n", ' ');
break;
case 22:
sprintf (out, "%6cRADF was appended successfully\n", ' ');
break;
default:
return 1;
break;
}
fputs (out, output);
return 0;
}

14
mm_radf/src/messages.h Normal file
View File

@ -0,0 +1,14 @@
/**
* @file
*/
#ifndef MESSAGES_H
#define MESSAGES_H
/**
* @fn message
*/
int message (const int, const int, const char *, FILE *);
#endif /* MESSAGES_H */

302
mm_radf/src/radf.c Normal file
View File

@ -0,0 +1,302 @@
/**
* @file
*/
#include <math.h>
// pi
#if !defined __USE_BSD && !defined __USE_XOPEN
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif /* __USE_BSD && __USE_XOPEN */
#endif /* M_PI */
/**
* @fn search_rdf
*/
int search_rdf (const int num_atoms, const int *type_atoms, const int *label_mol,
const float *coords, const double r_min, const double r_max,
const double r_step, int *radf)
/**
* @brief function that searchs molecule for rdf massive
* @code
* search_rdf (num_atoms, type_atoms, label_mol, coords, r_min, r_max, r_step, radf);
* @endcode
*
* @param num_atoms number of atoms
* @param type_atoms massive of atom types
* @param label_mol massive of numbers of molecule for atoms
* @param coords massive of coordinates
* @param r_min minimal radius
* @param r_max maximal radius
* @param r_step radius step
* @param radf not normed RDF
*
* @return 0 - exit without errors
*/
{
float r;
int i, j, k;
/* r radius
*/
for (i=0; i<num_atoms; i++)
for (j=0; j<8*num_atoms; j++)
if (((type_atoms[i] == 0) && (type_atoms[j] == 1)) ||
((type_atoms[i] == 1) && (type_atoms[j] == 0)))
{
r = sqrt (pow ((coords[3*i+0] - coords[3*j+0]), 2) +
pow ((coords[3*i+1] - coords[3*j+1]), 2) +
pow ((coords[3*i+2] - coords[3*j+2]), 2));
if ((r >= r_min) && (r <= r_max))
{
k = (r - r_min) / r_step;
radf[k]++;
}
}
return 0;
}
/**
* @fn search_rdf_centr
*/
int search_rdf_centr (const int num_atoms, const int *type_atoms, const int *label_mol,
const float *coords, const double r_min, const double r_max,
const double r_step, int *radf)
/**
* @brief function that searchs molecule for rdf massive by centered coordinates
* @code
* search_rdf_centr (num_atoms, type_atoms, label_mol, coords, r_min, r_max, r_step,
* radf);
* @endcode
*
* @param num_atoms number of atoms
* @param type_atoms massive of atom types
* @param label_mol massive of numbers of molecule for atoms
* @param coords massive of coordinates
* @param r_min minimal radius
* @param r_max maximal radius
* @param r_step radius step
* @param radf not normed RDF
*
* @return 0 - exit without errors
* @return 1 - error in set center (missing atoms)
*/
{
float r, cur_coords[2][3];
int coef, cur_at, i, j, k, l;
/* cur_coords centered coordinates
* r radius
* coef ==1 if j<num_atoms, ==7 if j>=num_atoms
* cur_at number of founded atoms
*/
for (i=0; i<num_atoms; i++)
if ((type_atoms[i] == 0) || (type_atoms[i] == 3))
for (j=0; j<8*num_atoms; j++)
if (((type_atoms[i] == 0) && (type_atoms[j] == 3)) ||
((type_atoms[i] == 3) && (type_atoms[j] == 0)))
{
// set center for i-molecule
for (l=0; l<3; l++)
cur_coords[0][l] = coords[3*i+l];
cur_at = 1;
for (k=1; k<6; k++)
if (label_mol[i+k] == label_mol[i])
if ((type_atoms[i+k] == type_atoms[i] + 1) ||
(type_atoms[i+k] == type_atoms[i] + 2))
{
cur_at++;
for (l=0; l<3; l++)
cur_coords[0][l] += coords[3*(i+k)+l];
}
if (cur_at == 3)
for (l=0; l<3; l++)
cur_coords[0][l] = cur_coords[0][l] / cur_at;
else
return 1;
// set center for j-molecule
for (l=0; l<3; l++)
cur_coords[1][l] = coords[3*j+l];
cur_at = 1;
if (j < num_atoms)
coef = 1;
else
coef = 7;
for (k=1; k<6; k++)
if (label_mol[j+coef*k] == label_mol[j])
if ((type_atoms[j+coef*k] == type_atoms[j] + 1) ||
(type_atoms[j+coef*k] == type_atoms[j] + 2))
{
cur_at++;
for (l=0; l<3; l++)
cur_coords[1][l] += coords[3*(j+coef*k)+l];
}
if (cur_at == 3)
for (l=0; l<3; l++)
cur_coords[1][l] = cur_coords[1][l] / cur_at;
else
return 1;
r = sqrt (pow ((cur_coords[0][0] - cur_coords[1][0]), 2) +
pow ((cur_coords[0][1] - cur_coords[1][1]), 2) +
pow ((cur_coords[0][2] - cur_coords[1][2]), 2));
if ((r >= r_min) && (r <= r_max))
{
k = (r - r_min) / r_step;
radf[k]++;
}
}
return 0;
}
/**
* @fn search_radf
*/
int search_radf (const int num_atoms, const int *type_atoms, const int *label_mol,
const float *coords, const double r_min, const double r_max,
const double r_step, const double ang_min, const double ang_max,
const double ang_step, int *radf)
/**
* @brief function that searchs molecule for radf massive
* @code
* search_radf (num_atoms, type_atoms, label_mol, coords, r_min, r_max, r_step,
* ang_min, ang_max, ang_step, radf);
* @endcode
*
* @param num_atoms number of atoms
* @param type_atoms massive of atom types
* @param label_mol massive of numbers of molecule for atoms
* @param coords massive of coordinates
* @param r_min minimal radius
* @param r_max maximal radius
* @param r_step radius step
* @param ang_min minimal angle
* @param ang_max maximal angle
* @param ang_step anlge step
* @param radf not normed RADF
*
* @return 0 - exit without errors
* @return 1 - error in set center (missing atoms)
*/
{
float ang, cos_ang, cur_coords[2][3], normal[2][3], r;
int coef, cur_at, i, j, k, l, n, n_max, plane_index[6];
/* ang angle
* cur_coords centered coordinates
* normal normal to planes
* r radius
* coef ==1 if j<num_atoms, ==7 if j>=num_atoms
* cur_at number of founded atoms
* n_max range of angles
* plane_index atoms forming plane
*/
for (i=0; i<num_atoms; i++)
if ((type_atoms[i] == 0) || (type_atoms[i] == 3))
for (j=0; j<8*num_atoms; j++)
if (((type_atoms[i] == 0) && (type_atoms[j] == 3)) ||
((type_atoms[i] == 3) && (type_atoms[j] == 0)))
{
// set center for i-molecule
plane_index[type_atoms[i]] = i;
for (l=0; l<3; l++)
cur_coords[0][l] = coords[3*i+l];
cur_at = 1;
for (k=1; k<6; k++)
if (label_mol[i+k] == label_mol[i])
if ((type_atoms[i+k] == type_atoms[i] + 1) ||
(type_atoms[i+k] == type_atoms[i] + 2))
{
cur_at++;
plane_index[type_atoms[i+k]] = i + k;
for (l=0; l<3; l++)
cur_coords[0][l] += coords[3*(i+k)+l];
}
if (cur_at == 3)
for (l=0; l<3; l++)
cur_coords[0][l] = cur_coords[0][l] / cur_at;
else
return 1;
// set center for j-molecule
plane_index[type_atoms[j]] = j;
for (l=0; l<3; l++)
cur_coords[1][l] = coords[3*j+l];
cur_at = 1;
if (j < num_atoms)
coef = 1;
else
coef = 7;
for (k=1; k<6; k++)
if (label_mol[j+coef*k] == label_mol[j])
if ((type_atoms[j+coef*k] == type_atoms[j] + 1) ||
(type_atoms[j+coef*k] == type_atoms[j] + 2))
{
cur_at++;
plane_index[type_atoms[j+coef*k]] = j + coef * k;
for (l=0; l<3; l++)
cur_coords[1][l] += coords[3*(j+coef*k)+l];
}
if (cur_at == 3)
for (l=0; l<3; l++)
cur_coords[1][l] = cur_coords[1][l] / cur_at;
else
return 1;
r = sqrt (pow ((cur_coords[0][0] - cur_coords[1][0]), 2) +
pow ((cur_coords[0][1] - cur_coords[1][1]), 2) +
pow ((cur_coords[0][2] - cur_coords[1][2]), 2));
// define planes
normal[0][0] = (coords[3*plane_index[1]+1] - coords[3*plane_index[0]+1]) *
(coords[3*plane_index[2]+2] - coords[3*plane_index[0]+2]) -
(coords[3*plane_index[1]+2] - coords[3*plane_index[0]+2]) *
(coords[3*plane_index[2]+1] - coords[3*plane_index[0]+1]);
normal[0][1] = (coords[3*plane_index[1]+2] - coords[3*plane_index[0]+2]) *
(coords[3*plane_index[2]+0] - coords[3*plane_index[0]+0]) -
(coords[3*plane_index[1]+0] - coords[3*plane_index[0]+0]) *
(coords[3*plane_index[2]+2] - coords[3*plane_index[0]+2]);
normal[0][2] = (coords[3*plane_index[1]+0] - coords[3*plane_index[0]+0]) *
(coords[3*plane_index[2]+1] - coords[3*plane_index[0]+1]) -
(coords[3*plane_index[1]+1] - coords[3*plane_index[0]+1]) *
(coords[3*plane_index[2]+0] - coords[3*plane_index[0]+0]);
normal[1][0] = (coords[3*plane_index[4]+1] - coords[3*plane_index[3]+1]) *
(coords[3*plane_index[5]+2] - coords[3*plane_index[3]+2]) -
(coords[3*plane_index[4]+2] - coords[3*plane_index[3]+2]) *
(coords[3*plane_index[5]+1] - coords[3*plane_index[3]+1]);
normal[1][1] = (coords[3*plane_index[4]+2] - coords[3*plane_index[3]+2]) *
(coords[3*plane_index[5]+0] - coords[3*plane_index[3]+0]) -
(coords[3*plane_index[4]+0] - coords[3*plane_index[3]+0]) *
(coords[3*plane_index[5]+2] - coords[3*plane_index[3]+2]);
normal[1][2] = (coords[3*plane_index[4]+0] - coords[3*plane_index[3]+0]) *
(coords[3*plane_index[5]+1] - coords[3*plane_index[3]+1]) -
(coords[3*plane_index[4]+1] - coords[3*plane_index[3]+1]) *
(coords[3*plane_index[5]+0] - coords[3*plane_index[3]+0]);
cos_ang = (normal[0][0] * normal[1][0] + normal[0][1] * normal[1][1] + normal[0][2] * normal[1][2]) /
(sqrt (pow ((normal[0][0]), 2) + pow ((normal[0][1]), 2) + pow ((normal[0][2]), 2)) *
sqrt (pow ((normal[1][0]), 2) + pow ((normal[1][1]), 2) + pow ((normal[1][2]), 2)));
ang = 180 / M_PI * acos (cos_ang);
if (ang > 90)
ang = 180 - ang;
n_max = (ang_max - ang_min) / ang_step;
if ((r >= r_min) && (r <= r_max))
if ((ang >= ang_min) && (ang <= ang_max))
{
k = (r - r_min) / r_step;
n = (ang - ang_min) / ang_step;
radf[n_max*k+n]++;
}
}
return 0;
}

26
mm_radf/src/radf.h Normal file
View File

@ -0,0 +1,26 @@
/**
* @file
*/
#ifndef RADF_H
#define RADF_H
/**
* @fn search_rdf
*/
/**
* @fn search_rdf_centr
*/
/**
* @fn search_radf
*/
int search_rdf (const int, const int *, const int *, const float *, const double,
const double, const double, int *);
int search_rdf_centr (const int, const int *, const int *, const float *, const double,
const double, const double, int *);
int search_radf (const int, const int *, const int *, const float *, const double,
const double, const double, const double, const double, const double,
int *);
#endif /* RADF_H */

137
mm_radf/src/radf_proc.c Normal file
View File

@ -0,0 +1,137 @@
/**
* @file
*/
#include <math.h>
#include <stdio.h>
// pi
#if !defined __USE_BSD && !defined __USE_XOPEN
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif /* __USE_BSD && __USE_XOPEN */
#endif /* M_PI */
/**
* @fn print_result
*/
int print_result (const char *output, const int matrix, const int mode,
const int step, const int num_atoms, const double r_min,
const double r_max, const double r_step, const double ang_min,
const double ang_max, const double ang_step, const float *cell,
const int *radf)
/**
* @brief function that print result to output file
* @code
* print_result (output, matrix, mode, step, num_atoms, r_min, r_max, r_step, ang_min,
* ang_max, ang_step, cell, radf);
* @endcode
*
* @param output output file name
* @param matrix status of matrix-mode
* @param mode 1 - if RDF, 2 - if RDF for center mass, 3 - if RADF
* @param step $(to - from + 1)
* @param num_atoms number of atoms
* @param r_min minimal radius
* @param r_max maximal radius
* @param r_step radius step
* @param ang_min minimal angle
* @param ang_max maximal angle
* @param ang_step angle step
* @param cell cell size
* @param radf not normed RADF
*
* @return 0 - exit without errors
*/
{
float ang, dv, norm, r, ro;
int i, j, k, n;
FILE *f_out;
/* ang angle
* dv volume
* norm normed RDF
* r radius
* ro density
* f_out output file
*/
f_out = fopen (output, "a");
fprintf (f_out, "SUMMARY STATISTIC\n");
switch (mode)
{
case 0:
if (matrix == 0)
fprintf (f_out, "| r | dV | RDF | RDFnorm |\n------------------------------------------\n");
k = (r_max - r_min) / r_step;
for (i=0; i<k; i++)
{
r = r_min + (0.5 + i) * r_step;
dv = 4 * M_PI * pow (r, 2) * r_step;
ro = num_atoms / (2 * cell[0] * cell[1] * cell[2]);
norm = radf[i] / (dv * ro * num_atoms * step);
if (matrix == 0)
fprintf (f_out, " %9.4f %10.4e %9i %9.6f \n", r, dv, radf[i] / 2, norm);
else
fprintf (f_out, " %9.4f %9.6f\n", r, norm);
}
fprintf (f_out, "------------------------------------------\n");
break;
case 1:
if (matrix == 0)
fprintf (f_out, "| r | dV | RDF | RDFnorm |\n------------------------------------------\n");
k = (r_max - r_min) / r_step;
for (i=0; i<k; i++)
{
r = r_min + (0.5 + i) * r_step;
dv = 4 * M_PI * pow (r, 2) * r_step;
ro = num_atoms / (3 * 2 * cell[0] * cell[1] * cell[2]);
norm = radf[i] / (dv * ro * num_atoms / 3 * step);
if (matrix == 0)
fprintf (f_out, " %9.4f %10.4e %9i %9.6f \n", r, dv, radf[i] / 2, norm);
else
fprintf (f_out, " %9.4f %9.6f\n", r, norm);
}
fprintf (f_out, "------------------------------------------\n");
break;
case 2:
if (matrix == 0)
fprintf (f_out, "| r | ang | dV | RDF | RDFnorm |\n----------------------------------------------------\n");
k = (r_max - r_min) / r_step;
n = (ang_max - ang_min) / ang_step;
if (matrix == 1)
{
fprintf (f_out, " r\\ang ");
for (j=0; j<n; j++)
fprintf (f_out, " %9.2f", ang_min+(0.5+j)*ang_step);
fprintf (f_out, "\n");
}
for (i=0; i<k; i++)
{
if (matrix == 1)
fprintf (f_out, " %9.4f", r);
for (j=0; j<n; j++)
{
r = r_min + (0.5 + i) * r_step;
ang = ang_min + (0.5 + j) * ang_step;
dv = 4 * M_PI * pow (r, 2) * sin (M_PI * ang / 180) * r_step * (M_PI * ang_step / 180);
ro = num_atoms / (3 * 2 * cell[0] * cell[1] * cell[2]);
norm = radf[n*i+j] / (dv * ro * num_atoms / 3 * step);
if (matrix == 0)
fprintf (f_out, " %9.4f %9.2f %10.4e %9i %9.6f \n", r, ang, dv, radf[n*i+j] / 2, norm);
else
fprintf (f_out, " %9.6f", norm);
}
if (matrix == 1)
fprintf (f_out, "\n");
}
fprintf (f_out, "----------------------------------------------------\n");
break;
}
fclose (f_out);
return 0;
}

16
mm_radf/src/radf_proc.h Normal file
View File

@ -0,0 +1,16 @@
/**
* @file
*/
#ifndef RADF_PROC_H
#define RADF_PROC_H
/**
* @fn print_result
*/
int print_result (const char *, const int, const int, const int, const int, const double,
const double, const double, const double, const double, const double,
const float *, const int *);
#endif /* RADF_PROC_H */