Commit 3723f56d authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Implement essential print_model_summary() function

parent 92e1b800
Loading
Loading
Loading
Loading
+60 −6
Original line number Diff line number Diff line
@@ -31,7 +31,7 @@ from sys import argv
## \cond
__version__ = "0.10.10"
int_reg = re.compile(r'[-+]?[0-9]+')
number_reg = re.compile(r'[-+]?[0-9]+\.[0-9]+(E[-+]?[0-9]+)?')
number_reg = re.compile(r'[-+]?[0-9]+\.[0-9]+([eE][-+]?[0-9]+)?')
## \endcond

## \brief Main execution code.
@@ -57,7 +57,7 @@ def main():
        else:
            model = scan_model(config['edfb_file'], config['geom_file'])
            if (model is not None):
                print(model)
                print_model_summary(model)
            else:
                result = 1
    return result
@@ -128,6 +128,58 @@ def print_help():
    print("--version             Print script version and exit.      ")
    print("                                                          ")

## \brief Print a summary of model properties.
#
#  This function provides a summary of useful information concerning the
#  radii of the particle monomers and of the equivalent mass sphere, to
#  assist in the selection of the proper starting orders.
#
#  \param[in] model: `dict` A model description dictionary.
def print_model_summary(model):
    avgX = 0.0
    avgY = 0.0
    avgZ = 0.0
    Rmin = 0.0
    Rmax = 0.0
    Reqm = 0.0
    R3tot = 0.0
    Rcirc = 0.0
    square_farthest = 0.0
    rad_farthest = 0.0
    for i in range(model['nsph']):
        sph_type_index = model['vec_types'][i] - 1
        ros = model['ros'][sph_type_index]
        avgX += model['vec_x'][i]
        avgY += model['vec_y'][i]
        avgZ += model['vec_z'][i]
        R3tot += math.pow(ros, 3.0)
        if (ros > Rmax):
            Rmax = ros
        if (Rmin == 0.0 or ros < Rmin):
            Rmin = ros
    Reqm = math.pow(R3tot, 1.0 / 3.0)
    avgX /= model['nsph']
    avgY /= model['nsph']
    avgZ /= model['nsph']
    for i in range(model['nsph']):
        sph_type_index = model['vec_types'][i] - 1
        ros = model['ros'][sph_type_index]
        dX = model['vec_x'][i] - avgX
        dY = model['vec_y'][i] - avgY
        dZ = model['vec_z'][i] - avgZ
        square_range = dX * dX + dY * dY + dZ * dZ + ros * ros
        if (square_range > square_farthest):
            square_farthest = square_range
            rad_farthest = ros
    Rcirc = math.sqrt(square_farthest) + rad_farthest
    if (Rmax == Rmin):
        print("INFO: monomer radius is Rsph = %.5em"%Rmin)
    else:
        print("INFO: smallest monomer radius is Rmin = %.5em"%Rmin)
        print("INFO: largest monomer radius is Rmax = %.5em"%Rmax)
    print("INFO: equivalent volume radius is Reqv = %.5em"%Reqm)
    print("INFO: minimum encircling radius is Rcirc = %.5em"%Rcirc)

## \brief Scan the calculation model.
#
#  The computational costs depennd the characteristics of the model. This
@@ -145,6 +197,7 @@ def scan_model(edfb_name, geom_name):
        'le': 0,
        'configurations': 0,
        'ros': [],
        'vec_types': [],
        'vec_x': [],
        'vec_y': [],
        'vec_z': []
@@ -216,6 +269,7 @@ def scan_model(edfb_name, geom_name):
            found_spheres += len(array_values)
            for ci in range(len(array_values)):
                type_id = int(array_values[ci].group())
                model['vec_types'].append(type_id)
                if (type_id > last_configuration):
                    last_configuration = type_id
                    configurations += 1
@@ -302,6 +356,7 @@ def scan_model(edfb_name, geom_name):
            file_line = edfb_file.readline()
            file_line = edfb_file.readline()
            read_lines += 2
            model['vec_types'].append(1)
            iter_values = int_reg.finditer(file_line)
            array_values = iter_to_array(iter_values)
            if (len(array_values) < 1):
@@ -373,6 +428,7 @@ def scan_model(edfb_name, geom_name):
                found_spheres += len(array_values)
                for ci in range(len(array_values)):
                    type_id = int(array_values[ci].group())
                    model['vec_types'].append(type_id)
                    if (type_id > last_configuration):
                        last_configuration = type_id
                        configurations += 1
@@ -420,13 +476,13 @@ def scan_model(edfb_name, geom_name):
        array_values = iter_to_array(iter_values)
        if (len(array_values) != 9):
            print("ERROR: %s is not a valid geometry configuration file!"%geom_name)
            print("   INVALID LINE: \"%s\""%file_line[:-1])
            print("   at line %d in %s."%(read_lines, geom_name))
            geom_file.close()
            return None
        nsph = int(array_values[0].group())
        if (nsph != model['nsph']):
            print("ERROR: %s is not consistent with %s!"%(geom_name, edfb_name))
            print("   INVALID LINE: \"%s\""%file_line[:-1])
            print("   at line %d in %s."%(read_lines, geom_name))
            geom_file.close()
            return None
        model['li'] = int(array_values[1].group())
@@ -439,8 +495,6 @@ def scan_model(edfb_name, geom_name):
            array_values = iter_to_array(iter_values)
            if (len(array_values) != 3):
                print("ERROR: %s is not consistent with %s!"%(geom_name, edfb_name))
                print("   INVALID LINE: \"%s\""%file_line[:-1])
                print("   at line %d in %s."%(read_lines, geom_name))
                geom_file.close()
                return None
            x = float(array_values[0].group())