#include "M_DynaMix_input.h"

#include "Log.h"

#include "Text.h"
#include "User.h"
#include "Path.h"
#include "M_DynaMix_input.h"
#include "Boundary_conditions.h"
#include "Model.h"
#include "Model_kit.h"
#include "Defs.h"

#include "M_DynaMix_mol_file.h" //fix remove with global_sigma_rule_


#include <fstream>
#include <sstream>
#include <string>

namespace MM
{
Duration_unit M_DynaMix_input::
time_step () const 
{
    return time_step_;
}

void M_DynaMix_input::
set_time_step (Duration_unit l) 
{
    time_step_ = l;
}


M_DynaMix_input::
M_DynaMix_input (Model & model)
:
    model_(&model),
//    file_name_                  ("M.DynaMix.input")
//    out_(new File_text_output   (file_name_)),
    source_file_name_           ("MDynaMix.input"),//fix
    output_control_             (5),//7
    //input_file_name_            ("md.input"),
    main_file_name_             ("MDynaMix"),
    //molecular_database_path_    (Path::MDynaMix_database().c_str()),//fix remove or move to
    read_from_restart_file_     (false),

    dump_restart_file_                                  (true),
    check_only_                                         (false),
    zero_counter_of_cpu_time_                           (true),

    constant_temperature_                               (true),
    constant_pressure_                                  (false),  //false
    anizotropic_NPT_                                    (false),
    
    temperature_                                        (Constant::normal_temperature),
    density_                                            (1.), //1.03 .6
    pressure_                                           (Constant::normal_pressure),
    //box_x_                                              (0.),
    //box_y_                                              (0.),
    //box_z_                                              (0.),
    cell_type_                                          (0),
    //boundary_conditions_                        (Boundary_conditions::vacuum),
    //box_type_                                   (Box_type::automatic),
    
    time_step_                                          (0.002), //(2.0),

    small_steps_in_one_long_                            (10),    //10
//    total_long_steps_                                   (1000),
//    total_long_steps_                                   (10),
    steps_for_intermediate_averaging_                   (1000),
    take_averages_each___steps_                         (10), // SRP ?
    dump_restart_file_after___steps_                    (500),
    
    t_termostat_                                        (30.),
    p_termostat_                                        (700.),
    is_simple_velosity_scaling_                         (false),
    delta_t_                                            (20.),
    
    r_cutoff_                                           (12.), //10
    r_cutoff_fast_force_                                (5.),
    check_neigbours_after___step_                       (10),
    
//    alpha_r_                                            (3.14159256),
//    fexp_                                               (9.81),
    alpha_r_                                            (2.8),//2.7
    fexp_                                               (9.), //7.
    
    recalculate_list_of_intramolecular_interactions_    (true),
    
    constrain_dynamics_                                 (false),
    //constrain_dynamics_                                 (true),
    tolerance_parameter_                                (1.e-4),
    //tolerance_parameter_                                (1.e-2),

    initial_state_                                      (0),
    set_velocities_to_0_                                (false),//false
    hole_radius_                                        (10.),
    keep_initial_configuration_                         (false),
    permissible_deviation_                              (3.),
    file_name_                                          ("fixed.atoms"),

    change_temperature_at_restart_                      (false),
    change_density_at_restart_                          (false),

    final_averaging_after___intermediate_averaging_     (1), //6
    dump_XMOL_config_file_                              (true),

    //cut_large_forces_                                   (false),
    cut_large_forces_                                   (true),
    max_level_of_large_forces_                          (0.1),
    //max_level_of_large_forces_                          (5.e-5),//1.e-4

    calculate_RDF_                                      (false),//true
    dump_RDF_                                           (false),//true
    read_restart_RDF_file_                              (false),
    is_RDF_for_all_sites_                               (false),
    cutoff_RDF_                                         (10.),
    RDF_resolution_                                     (200),

    dump_trajectory_format_                             (2),
    dump_trajectory_interval_                           (1.e-14),
    number_of_config_in_a_trajectory_                   (500),

    calculate_tcf_                                      (false),
    restart_tcf_                                        (false),
    dump_tsf_                                           (false),
    NSTEG_                                              (200),
    JUMP_                                               (5),

    line_counter_(0)
{
    Duration_unit duration;
    duration.set_step_length_ps (0.002);
    duration.set_ps             (100.);
    set_time_step (duration);

    //model_.kit().boundary_conditions().set_type ("Box");
    //model_.kit().boundary_conditions().set_type ("Vacuum");
    //Model_kit & kit = model_.kit();

    M_DynaMix_molecular_type H2O ("H2O");
    H2O.number_of_molecules_ = 246;
    H2O.intramolecular_potential_type_ = 2;
    molecular_types_.push_back (H2O);
    
    M_DynaMix_molecular_type Na ("Na+");
    Na.number_of_molecules_ = 5;
    molecular_types_.push_back (Na);
    
    M_DynaMix_molecular_type Cl ("Cl-");
    Cl.number_of_molecules_ = 5;
    molecular_types_.push_back (Cl);
}

void M_DynaMix_input::
setup_dependency ()
{
    /*if (box_type_ == manual)
    {
    }
    else if (box_type_ == from_crystal)
    {
    }
    else if (box_type_ == automatic)
    {
    }
    else FLAW ("No such box type.");//*/

    /*Boundary_conditions const & boundary_conditions 
        = model_.kit().boundary_conditions();

    if (boundary_conditions.type_name() == "Vacuum")
    {
        box_x_ = 0;
        box_y_ = 0;
        box_z_ = 0;
    }
    else if (boundary_conditions_ == box)
    {
    }
    else FLAW ("No such boundary conditions.");//*/
}

void M_DynaMix_input::
set_vacuum ()
{
    //Boundary_conditions & boundary_conditions 
    //    = model_.kit().boundary_conditions();
    Model_kit & kit                           = model().kit();
    Boundary_conditions & boundary_conditions = kit.boundary_conditions();

    //boundary_conditions.set_type ("Box");
    //boundary_conditions.box ().set (0,0,0, 1,1,1);
    boundary_conditions.set_type ("Vacuum");
    //option_->box_x_   = 0;
    //option_->box_y_   = 0;
    //option_->box_z_   = 0;
//    density_ = 0;
    //alpha_r_ = 0;
    //fexp_    = 0;
}

void M_DynaMix_input::
set_box ()
{
    Boundary_conditions & boundary_conditions 
        = model().kit().boundary_conditions();

    boundary_conditions.set_type ("Box");
    //option_->box_x_   = 10;//fix
    //option_->box_y_   = 10;
    //option_->box_z_   = 10;
//    density_ = 1;
}

double M_DynaMix_input::
box_x ()
{
    REQUIRE ("Box exist", 
        model().kit().boundary_conditions().type_name() == "Box");

    return model().kit().boundary_conditions().box().a_mod_rectangle();
}

double M_DynaMix_input::
box_y ()
{
    REQUIRE ("Box exist", 
        model().kit().boundary_conditions().type_name() == "Box");

    return model().kit().boundary_conditions().box().b_mod_rectangle();
}

double M_DynaMix_input::
box_z ()
{
    REQUIRE ("Box exist", 
        model().kit().boundary_conditions().type_name() == "Box");

    return model().kit().boundary_conditions().box().c_mod_rectangle();
}

void M_DynaMix_input::
prepare_input (Text_output & input)
{
    to_user(). status ("Prepare input file");

    setup_dependency ();

    write_output_control                                    (input);
    write_output_file_name                                  (input);
    write_molecular_database_path                           (input);
    write_read_from_restart_file                            (input);
    write_dump_restart_file                                 (input);
    write_check_only                                        (input);
    write_zero_counter_of_cpu_time                          (input);
    write_type_of_statistical_ensemble                      (input);
    write_number_of_molecular_types                         (input);
    write_molecular_types                                   (input);
    write_number_of_molecules                               (input);
    write_nonbonded_intramolecular_interactions             (input);
    write_interactions_1_4                                  (input);
    write_scaling_factors_for_1_4_LJ                        (input);
    write_scaling_factors_for_1_4_electrostatics            (input);
    write_intramolecular_potential_type                     (input);
    write_rules_for_initial_box_size_or_density             (input);
    write_time_steps                                        (input);
    write_nose_termostat_parameters                         (input);
    write_cut_offs                                          (input);
    write_ewald_parameters                                  (input);
    write_which_of_the_molecules_move                       (input);
    write_recalculate_list_of_intramolecular_interactions   (input);
    write_constrain_dynamics                                (input);
    write_initial_state                                     (input);
    write_initial_states_for_atom_group_types               (input);
    write_initial_state_configuration                       (input);
    write_change_temparature_or_density_after_restart       (input);
    write_final_averaging                                   (input);
    write_cut_large_forces                                  (input);
    write_RDFs                                              (input);
    write_dump_trajectory                                   (input);
    write_TCF_calculations                                  (input);
    write_optional_parameters                               (input);
    write_RDFs_additional                                   (input);

//    input.line ("stub12345###");

    input.line ("");
    to_user(). status ("Input file prepared");
}

void M_DynaMix_input::
prepare_input_5 (Text_output & out)
{
    int i,j;

    to_user(). status ("Prepare input file");
    out.line ("# Input file for MDynaMxi v.>=5.0 generated by m.DynaMix.shell.");
    out.line ("#______________________________________________________________");
    out.line ("#");
    out.line ("# Basic setup");
    out.line ("#");

    out.text ("Main_filename          ");
    out.line (main_file_name_.c_str());
    //out.line ("md5");

    out.text ("Verbose_level          ");
    out.line (output_control_);

    out.text ("Path_DB                ");
    out.line (Path::MDynaMix () + "moldb");
    //out.line (""); //fix UNIX_path

    out.line ("Visual                 yes"); 

    out.line ("#");
    out.line ("# Restart control");
    out.line ("#");

    out.text ("Check_only             ");
    if (check_only_)                    out.line ("yes");
    else                                out.line ("no");

    out.text ("Read_restart           ");
    if (read_from_restart_file_)    out.line ("yes");
    else                            out.line ("no");

    if (dump_restart_file_)
    {
        out.text ("Dump_restart           ");
        out.line (dump_restart_file_after___steps_);
    }
    
    out.text ("Change_T               ");
    if (change_temperature_at_restart_) out.line ("yes");
    else                                out.line ("no");

    out.text ("Change_V               ");
    if (change_density_at_restart_)     out.line ("yes");
    else                                out.line ("no");

    //out.text ("Zero_CPU            ");
    //if (zero_counter_of_cpu_time_)      out.line ("yes");
    //else                                out.line ("no");

    out.text ("Zero_vel               ");
    if (set_velocities_to_0_)           out.line ("yes");
    else                                out.line ("no");

    out.line ("#");
    out.line ("# System");
    out.line ("#");

    out.text ("Molecule_types         ");
    out.line (molecular_types_.size());

    for (i=0;  i<molecular_types_.size();  ++i)
    {
        out.text (molecular_types_[i].name_);
        
        for (j=molecular_types_[i].name_.length();  j<22;  ++j)
            out.text (" ");
        out.text (" ");

        out.line (molecular_types_[i].number_of_molecules_);
    }
    
    if (model().kit().boundary_conditions().type_name() != "Vacuum")
    {
        out.text ("PBC                    ");
        out.line ("rect");
    }

    if (model().kit().boundary_conditions().type_name() == "Box") //fix
    {
        out.text ("Box                    ");
        out.text (box_x());
        out.text ("  ");
        out.text (box_y());
        out.text ("  ");
        out.line (box_z());
    }
    else 
    {
        out.line ("Density                0");
    }

    out.line ("#");
    out.line ("# Ensemble");
    out.line ("#");

    if (is_simple_velosity_scaling_)
    {
        out.text ("Velocity_scaling       ");
        out.text (temperature_);
        out.text (" ");
        out.line (delta_t_);
    }
    else
    {
        out.text ("Nose_thermostat        ");
        out.text (temperature_);
        out.text (" ");
        out.line (t_termostat_);
    }

    out.line ("Separate_thermostating yes");        //fix make control
    out.line ("COM_check              yes 100");    //fix make control

    if (constant_pressure_)
    {
        out.text ("Barostate_NH           ");
        out.text (pressure_);
        out.text ("   ");
        out.line (p_termostat_);
    }
    if (anizotropic_NPT_)           
        out.line ("Barostate_anisotropic  yes");

    out.line ("#");
    out.line ("# Molecular Dynamics");
    out.line ("#");

    out.text ("Time_step              ");   
    out.line (time_step_.step_length_ps() * 1000);

    out.text ("Number_steps           ");   
    //out.line (total_long_steps_);
    out.line ((int)time_step_.steps());

    if (constrain_dynamics_)
    {
        out.text ("Constrain              ");   
        out.text (tolerance_parameter_);
        out.text ("  "); 
        for (i=0;  i<molecular_types_.size();  ++i)
        {
            out.text (" "); 
            out.text ("1");
        }
        out.line ("");
    }
    else
    {
        out.text ("Double_timestep        ");   
        out.line (small_steps_in_one_long_);
    }

    out.text ("R_cutoff               ");   //fix auto
    out.line (r_cutoff_);           
    /*if (model_.kit().boundary_conditions().type_name() == "Box" && 
        model_.kit().boundary_conditions().is_auto_box ())
    {
        double L = (box_x()+box_y()+box_z())/3.; //fix
        if (L < 20 || L > 60)
            out.text (r_cutoff_);
        else
            out.text (0.1*(L-10.) + 10.);
    }
    else
        out.text (r_cutoff_);//*/

    out.text ("R_short                ");   
    out.line (r_cutoff_fast_force_);

    out.text ("Neighbour_list         ");   
    out.line (check_neigbours_after___step_);

    out.text ("Electrostatics         "); 
    if (model().kit().boundary_conditions().type_name() == "Vacuum")
    {
        out.line ("Cutoff");
    }
    else
    {
        out.text ("Ewald   ");
        out.text (alpha_r_);
        out.text (" ");
        out.line (fexp_);
    }

    if (cut_large_forces_)
    {
        out.text ("Cut_forces             ");
        out.line (max_level_of_large_forces_);
    }

    if (global_sigma_rule() == "geometrical")
    {
        out.line ("Combination_rule        Sigma_geom");
    }

    out.line ("#");
    out.line ("# Startup");
    out.line ("#");

    if (!read_from_restart_file_)
    {
        out.line ("Startup                xmol");
    }
     
    out.line ("Gather");

    out.line ("#");
    out.line ("# Properties");
    out.line ("#");

    out.text ("Output                 ");
    out.line (take_averages_each___steps_);
    
    out.text ("Serie_average          ");
    out.line (steps_for_intermediate_averaging_);
    
    out.text ("Average_from           ");
    out.line (final_averaging_after___intermediate_averaging_);
    
    if (dump_XMOL_config_file_)
        out.line ("Dump_XMOL              yes");
    else
        out.line ("Dump_XMOL              no");
    
    if (dump_trajectory_format_ != 0)
    {
        out.text ("Trajectory             ");
        switch (dump_trajectory_format_)
        {
        case 1:     out.text ("bincrd  ");    break;
        case 2:     out.text ("asccrd  ");    break;
        case 3:     out.text ("binvel  ");    break; //fix
        case 4:     out.text ("ascvel  ");    break; //fix
        default:
            FLAW (Text("No dump trajectory format of ") 
                + dump_trajectory_format_ + "type.");
        }

        out.text (dump_trajectory_interval_ * 1e15);
        out.text ("  ");

        out.text (number_of_config_in_a_trajectory_);
        out.text ("  ");

        out.line ("all");
    }

    /////////////////////////////////////////////
    out.line ("end");
    out.line ("");
    to_user(). status ("Input file prepared");
}

void M_DynaMix_input::
write_output_control (Text_output & out)
{
    out.line ("#   Sample input file for the MD program, v.>=4.3");
    out.line ("#");
    out.line ("#   Simulated system consists of 246 flexible SPC water molecules,");
    out.line ("#   and 5 Na+Cl- ion pairs");
    out.line ("#");
    out.line ("#   Lines beginning with \"#\" are commentaries");
    out.line ("#   In this file the commentaries go before the corresponding parameters.");
    out.line ("#");
    out.line ("#   This file can be directly used as input (see directory sample for");
    out.line ("#   short version of the input file)");
    out.line ("#");
    out.line ("#");
    out.line ("#   Output control parameter");
    out.line ("#   Suitable valus 2-10. The less number, the less you see in the output");
    out.line ("#   parameters higher 7 used mostly for debug purposes  ");
    
    out.text ("  ");
    out.line (output_control_);
}

void M_DynaMix_input::
write_output_file_name (Text_output & out)
{
    out.line ("#");
    out.line ("#   Base file name for output files:");
    out.line ("#   Other files requested or created by the program have this name with");
    out.line ("#   various extensions");
    
    out.text ("  ");
    out.line (main_file_name_.c_str());
}

void M_DynaMix_input::
write_molecular_database_path (Text_output & out)
{
    out.line ("#");
    out.line ("#   Path to the molecular database");
    out.line ("#   This directory contains *.mol files wich describe the molecular");
    out.line ("#   structure and the force field");
    
    out.text ("  ");
    //out.line (molecular_database_path_.c_str());
    out.line ("./moldb/"); //fix UNIX_path
}

void M_DynaMix_input::
write_read_from_restart_file (Text_output & out)
{
    out.line ("#");                  
    out.line ("# The program creates and updates periodically a restart file which");
    out.line ("# contains configuration of the system and calculated averages.");
    out.line ("# The program can be interrupted and then continued from the restart file");
    out.line ("# without loosing any information.");
    out.line ("#  If \"Check only\" parameter is true, the program does not run");
    out.line ("# simulation. If it is a new run, the program only checks the input.");
    out.line ("# If it is continuation of the old run, the program gives calculated results.");
    out.line ("#");
    out.line ("#  Read from        Dump               Check          Zero counter");
    out.line ("# restart file?  restart file?         only?          of cpu time?");

    out.text ("  ");
    if (read_from_restart_file_)    out.text (".t.");
    else                            out.text (".f.");
}

void M_DynaMix_input::
write_dump_restart_file (Text_output & out)
{
    out.text ("                ");
    if (dump_restart_file_)     out.text (".t.");
    else                        out.text (".f.");
}

void M_DynaMix_input::
write_check_only (Text_output & out)
{
    out.text ("                ");
    if (check_only_)            out.text (".t.");
    else                        out.text (".f.");
}

void M_DynaMix_input::
write_zero_counter_of_cpu_time (Text_output & out)
{
    out.text ("                ");
    if (zero_counter_of_cpu_time_)  out.line (".t.");
    else                            out.line (".f.");
}

void M_DynaMix_input::
write_type_of_statistical_ensemble (Text_output & out)
{
    out.line ("#");
    out.line ("#  The type of statistical ensemble. \"Anizotropic NPT\" means separate");
    out.line ("#  pressure/volume controll in each direction. Use this option only if the");
    out.line ("#  system is really anizotropic (a piece of DNA, membrane, liquid crystall,"); 
    out.line ("#  etc). Note, that for \"constant pressure\"=.t. \"constant temperature\"");
    out.line ("#  must be also .true. (who knows what is NVP ensemble?).");
    out.line ("#  Constant temperature?   Constant pressure?    Anizotropic NPT?");

    out.text ("  ");
    if (constant_temperature_)      out.text (".t.");
    else                            out.text (".f.");

    out.text ("                      ");
    if (constant_pressure_)         out.text (".t.");
    else                            out.text (".f.");

    out.text ("                    ");
    if (anizotropic_NPT_)           out.line (".t.");
    else                            out.line (".f.");
}

void M_DynaMix_input::
write_number_of_molecular_types (Text_output & out)
{
    out.line ("#");
    out.line ("#   Number of molecular types:");

    out.text ("  ");
    out.line (molecular_types_.size());
}

void M_DynaMix_input::
write_molecular_types (Text_output & out)
{
    out.line ("#");
    out.line ("#  Types of molecules (models):");
    out.line ("#  The database directory (specified above) should contain files");
    out.line ("#  H2O.mol, Na+.mol,  Cl-.mol");
    out.line ("#  See README file about format of .mol files");

    for (int i=0;  i<molecular_types_.size();  ++i)
    {
        out.text ("  ");
        out.text (molecular_types_[i].name_);
    }
    out.line ("");
}

void M_DynaMix_input::
erase_molecular_types ()
{
    molecular_types_.clear ();
}

void M_DynaMix_input::
add_molecular_type (M_DynaMix_molecular_type & new_type)
{
    molecular_types_.push_back (new_type);
}

void M_DynaMix_input::
increment_molecular_type (const Text & type)
{
    for (int i=0;  i<molecular_types_.size();  ++i)
        if (molecular_types_[i].name_ == type)
    {
        ++molecular_types_[i].number_of_molecules_;
        return;
    }

    FLAW ("No such molecular type");
}

void M_DynaMix_input::
write_number_of_molecules (Text_output & out)
{
    out.line ("#");
    out.line ("#   For each molecular type:");
    out.line ("#   Number of molecules:");

    for (int i=0;  i<molecular_types_.size();  ++i)
    {
        out.text ("  ");
        out.text (molecular_types_[i].number_of_molecules_);
    }
    out.line ("");
}

void M_DynaMix_input::
write_nonbonded_intramolecular_interactions (Text_output & out)
{
    out.line ("#");
    out.line ("#    non-bonded intramolecular interactions:");
    out.line ("#   Calculate intramolecular potential (LJ and electrostatic)");
    out.line ("#   for non-bound atoms, separated by");
    out.line ("#   more than 3 covalent bonds. Normally, it should be .true.");
    out.line ("#   unimportant for small molecules)");

    for (int i=0;  i<molecular_types_.size();  ++i)
    {
        out.text ("  ");
        if (molecular_types_[i].nonbonded_intramolecular_interactions_)
                out.text (".t.");
        else    out.text (".f.");
    }
    out.line ("");
}

void M_DynaMix_input::
write_interactions_1_4 (Text_output & out)
{
    out.line ("#");
    out.line ("#    1 - 4 intramolecular interactions:");
    out.line ("#   Calculate LJ and electrostatic terms for 1-4 (i.e. separated by");
    out.line ("#   3 covalent bonds) intramolecular interactions");

    for (int i=0;  i<molecular_types_.size();  ++i)
    {
        out.text ("  ");
        if (molecular_types_[i].interactions_1_4_)  out.text (".t.");
        else                                        out.text (".f.");
    }
    out.line ("");
}

void M_DynaMix_input::
write_scaling_factors_for_1_4_LJ (Text_output & out)
{
    out.line ("#   Scaling factors for 1-4 LJ interactions");
    out.line ("#   These are 0 in AMBER, 1 in CHARMM and 0.25 in GROMOS"); //fix
    out.line ("#   Not important for small molecules");

    for (int i=0;  i<molecular_types_.size();  ++i)
    {
        out.text ("  ");
        out.text (molecular_types_[i].scaling_factors_for_1_4_LJ_);
    }
    out.line ("");
}

void M_DynaMix_input::
write_scaling_factors_for_1_4_electrostatics (Text_output & out)
{
    out.line ("#   Scaling factor for 1-4 electrostatic interactions");

    for (int i=0;  i<molecular_types_.size();  ++i)
    {
        out.text ("  ");
        out.text (molecular_types_[i].scaling_factors_for_1_4_electrostatics_);
    }
    out.line ("");
}

void M_DynaMix_input::
write_intramolecular_potential_type (Text_output & out)
{
    out.line ("#");
    out.line ("#    intramolecular potential type:");
    out.line ("#    Now it should be 0 for all molecules except the water");
    out.line ("#    For water 1 is \"harmonic\" and 2 \"anharmonic\" flexible SPC water");
    out.line ("#    (It says the program to use a special subroutine)");

    for (int i=0;  i<molecular_types_.size();  ++i)
    {
        out.text ("  ");
        out.text (molecular_types_[i].intramolecular_potential_type_);
    }
    out.line ("");
}

void M_DynaMix_input::
write_rules_for_initial_box_size_or_density (Text_output & out)
{
    out.line ("#");
    out.line ("#   Rules for initial box size / density:");
    out.line ("#    - If one of the box sizes is zero, the actual box size (cubic box)");
    out.line ("#   is defined by the density."); //fix how?
    out.line ("#    - If the density is also zero, the program run \"vacuum");
    out.line ("#   simulations\". Set Ewald parameters (below) to 0 in the case of");
    out.line ("#   vacuum simulations.");
    out.line ("#    - If all the three box sizes are not zero, they define initial box size");
    out.line ("#   and shape, and so the actual density");
    out.line ("#");
    out.line ("#   Cell type:");
    out.line ("#   0 - rectangular");
    out.line ("#   1 - truncated octahedron: cube of side BOXL centred in 0,0,0 with");
    out.line ("#           truncated corners:      |x|+|y|+|z| < 0.75*BOXL");
    out.line ("#   2 - hexagonal along Z axis");
    out.line ("#");
    out.line ("#   temperature(K)     density(g/cm**3)       pressure (atm)");

    out.text ("  ");
    out.text (temperature_);

    if (model().kit().boundary_conditions().type_name() == "Vacuum")
    {
        out.text ("                   0");
    }
    else
    {
        out.text ("                   ");
        out.text (density_);
    }

    out.text ("                   ");
    out.line (pressure_);


    out.line ("#        box size (A)                         cell type");
        out.text ("         ");

    if (model().kit().boundary_conditions().type_name() == "Box")
    {
        out.text (box_x());
        out.text ("  ");
        out.text (box_y());
        out.text ("  ");
        out.text (box_z());
    }
    else 
    {
        out.text ("0  0  0");
    }

    out.text ("                              ");
    out.line (cell_type_);
}

void M_DynaMix_input::
write_time_steps (Text_output & out)
{
    out.line ("#");
    out.line ("#  This is the long time step");
    out.line ("#  Time step (s)      Small steps in one long");

    out.text ("  ");
    //out.text (time_step_ / 1e15);
    out.text (time_step_.fs());

    out.text ("                   ");
    out.line (small_steps_in_one_long_);


    out.line ("#");
    out.line ("# Total (long)   Steps for interme-      take averages     dump restart file");
    out.line ("#  MD-steps      diate averaging         each .. steps      after .. steps");
    
    out.text ("  ");
    //out.text (total_long_steps_);
    out.text (time_step_.steps());
    
    out.text ("                    ");
    out.text (steps_for_intermediate_averaging_);
    
    out.text ("                    ");
    out.text (take_averages_each___steps_);
    
    out.text ("                    ");
    out.line (dump_restart_file_after___steps_);
}

void M_DynaMix_input::
write_nose_termostat_parameters (Text_output & out)
{
    out.line ("#");
    out.line ("# Nose termostat parameters:                 Meaningful in constant-energy");
    out.line ("#                                            simulations (if const.temp.=.f.):");
    out.line ("#  T-termostat param(fs)  P-termostat(fs)    Simple velosity   delta T (K)");
    out.line ("#                                             scaling?");

    out.text ("  ");
    out.text (t_termostat_);

    out.text ("                    ");
    out.text (p_termostat_);

    out.text ("                    ");
    if (is_simple_velosity_scaling_)    out.text (".t.");
    else                                out.text (".f.");

    out.text ("                    ");
    out.line (delta_t_);
}

void M_DynaMix_input::
write_cut_offs (Text_output & out) //fix GUI
{
    out.line ("#");
    out.line ("#  Rcutoff sets cut off radius for LJ and Real-space electrostatic forces");
    out.line ("#  Interaction inside Rcut-fast are recalculated each short time step,");
    out.line ("#  interactions between Rcut-fast and Rcutoff each long time step");
    out.line ("#");
    out.line ("#  Rcutoff(A)      Rcut-fast forces      check neigbours after .. steps");
    
    out.text ("  ");

    if (model().kit().boundary_conditions().type_name() == "Box" && 
        model().kit().boundary_conditions().is_auto_box ())
    {
        double L = (box_x()+box_y()+box_z())/3.; //fix
        if (L < 20 || L > 60)
            out.text (r_cutoff_);
        else
            out.text (0.1*(L-10.) + 10.);
    }
    else
        out.text (r_cutoff_);

    out.text ("                  ");
    out.text (r_cutoff_fast_force_);

    out.text ("                   ");
    out.line (check_neigbours_after___step_);
}

void M_DynaMix_input::
write_ewald_parameters (Text_output & out)
{
    out.line ("#===========================================");
    out.line ("#   Treatment of electrostatic interactions:");
    
    out.line ("#");
    out.line ("#   Ewald parameters:");
    out.line ("#   alpha and fexp are set from the conditions:");
    out.line ("#  erfc(alpha*R)=required precision of the real-space Ewald");
    out.line ("#  exp(-fexp)=required precision of reciprocal-space Ewald");
    out.line ("#    the rule of thumb:");
    out.line ("#    alpha*R                 fexp (m/s**2)");

    if (model().kit().boundary_conditions().type_name() == "Vacuum") 
    {
        out.line ("       0                       0");
    }
    else
    {
        out.text ("       ");
        out.text (alpha_r_);
 
        out.text ("                    ");
        out.line (fexp_);
    }

    out.line ("#   If alpha*R above is negative, the reaction field method is used with");
    out.line ("#   -alpha*R as dielectric permitivity, fexp is Debay screening length in Å");
    out.line ("#   (setting Debay length to 0 means infinite Debue length, i.e non-conducting");
    out.line ("#   solution)");
    out.line ("#   If alpha is zero (exactly: between -1 and 0), no special treatment of");
    out.line ("#   electrostatic interactions (simple spherical cutoff)");
    out.line ("#");
    out.line ("#   Be careful by playing with these parameters and understand what are you");
    out.line ("#   doing. Unproper setting can result in funny behavior of the system or too");
    out.line ("#   long computation time.");
}

void M_DynaMix_input::
write_which_of_the_molecules_move (Text_output & out)
{
    out.line ("#");
    out.line ("#===============================================");
    out.line ("#   which of the molecules move:");
    out.line ("#  if .f., the molecules are fixed, but interact with other molecules.");

    for (int i=0;  i<molecular_types_.size();  ++i)
    {
        out.text ("  ");
        if (molecular_types_[i].mobile_)    out.text (".t.");
        else                                out.text (".f.");
    }
    out.line ("");
}

//fix ???
void M_DynaMix_input::
write_recalculate_list_of_intramolecular_interactions (Text_output & out)
{
    out.line ("#");
    out.line ("#   recalculate the list of intramolecular interactions?");
    out.line ("#   can be set to .f., if you have a large molecule with"); 
    out.line ("#   a stable conformation");

    out.text ("  ");
    if (recalculate_list_of_intramolecular_interactions_)   out.line (".t.");
    else                                                    out.line (".f.");
}

void M_DynaMix_input::
write_constrain_dynamics (Text_output & out)
{
    out.line ("#");
    out.line ("#  If true, constrain dynamics with SHAKE algorithm for specified");
    out.line ("#  molecular specii will be used. It will keep all the bond lengths constant.");
    out.line ("#  No double time step alogorithm in this case");
    out.line ("#                                                 which molecules considered as rigid (1/0)");
    out.line ("#  Constrain dynamics?    tolerance parameter   (does not matter if Constrained dynamics = .f.)");

    out.text ("  ");
    if (constrain_dynamics_)    out.text (".t.");
    else                        out.text (".f.");
    
    out.text ("                     ");
    out.text (tolerance_parameter_);
    
    //fix
    out.text ("                ");
    out.line ("1 1 1 1 1 1 1 1 1");
}

void M_DynaMix_input::
write_initial_state (Text_output & out)
{
    out.line ("#");
    out.line ("#  Initial state (from -1 to 4)");
    out.line ("#  Values:");
    out.line ("#  -1  take initial center-of-mass coordinates of molecules from *.inp file");
    out.line ("#  0  take initial atom coordinates from *.inp file");
    out.line ("# /*    .inp file should be written in free  x y z format (one atom per line)");
    out.line ("#       The order of atoms and molecules:");
    out.line ("#       molecular type1                        molecular type 2");
    out.line ("#       mol1    mol2    mol3                    mol1    mol2    mol3");
    out.line ("#       at1 at2 at1 at2 at1 at2                at1 at2 at1 at2    ...   */");
    out.line ("#  1  start from FCC lattice");
    out.line ("#  2  set a cylindrical hole along z-axis. Put molecules with parameter 1");
    out.line ("#     specified on the line below into the hole according initial");
    out.line ("#     coordinates in .mol file, and distribute others out the hole");
    out.line ("#  3  the same for spherical hole");
    out.line ("#  4  start from cubic lattice");
    out.line ("#                              |             set velocities to 0?");
    out.line ("#                              |             (otherwise, Maxwel distribution,");
    out.line ("#                              |             or as it is at restart)");

    out.text ("  ");
    out.text (initial_state_);
    
    out.text ("                                                ");
    if (set_velocities_to_0_)   out.line (".t.");
    else                        out.line (".f.");
}

void M_DynaMix_input::
write_initial_states_for_atom_group_types (Text_output & out)
{
    out.line ("#");
    out.line ("#  Parameters for above in cases 2,3 (specify 0 or 1)");
    
    for (int i=0;  i<molecular_types_.size();  ++i)
    {
        out.text ("  ");
        out.text (molecular_types_[i].initial_state_);
    }
    out.line ("");
}

void M_DynaMix_input::
write_initial_state_configuration (Text_output & out)
{
    out.line ("#");
    out.line ("#                              If this parameter .t., the atoms described in");
    out.line ("#                              file \"fixed.atoms\" will be");
    out.line ("#                              put in a harmonic potential of this radius:");
    out.line ("#  Radius of a hole            Keep initial          Permissible     File");
    out.line ("#  for large molecules         configuration?        deviation?      name");

    out.text ("  ");
    out.text (hole_radius_);
    
    out.text ("                           ");
    if (keep_initial_configuration_)    out.text (".t.");
    else                                out.text (".f.");

    out.text ("                   ");
    out.text (permissible_deviation_);
    
    out.text ("             ");
    out.line (file_name_);
}

void M_DynaMix_input::
write_change_temparature_or_density_after_restart (Text_output & out)
{
    out.line ("#");
    out.line ("#  If you want to change temparature or density after restart.");
    out.line ("#  Change temperature at restart?      Change density?");

    out.text ("  ");
    if (change_temperature_at_restart_) out.text (".t.");
    else                                out.text (".f.");

    out.text ("                                  ");
    if (change_density_at_restart_) out.line (".t.");
    else                            out.line (".f.");
}

void M_DynaMix_input::
write_final_averaging (Text_output & out)
{
    out.line ("#");
    out.line ("#  When to start final averaging");
    out.line ("#  Final averaging after ..           Dump");
    out.line ("#  intermediate averaging            XMOL config. file?");

    out.text ("  ");
    out.text (final_averaging_after___intermediate_averaging_);

    out.text ("                                    ");
    if (dump_XMOL_config_file_)     out.line (".t.");
    else                            out.line (".f.");
}

void M_DynaMix_input::
write_cut_large_forces (Text_output & out)
{
    out.line ("#");
    out.line ("#  This have a sence if you have to start from a very bad configuration");
    out.line ("#  If total force on an atom exceed some level defined by the given");
    out.line ("#  parameter, the force will be cut to this level in the internal units");
    out.line ("#  Cut large forces?         Parameter");

    out.text ("  ");
    if (cut_large_forces_)      out.text (".t.");
    else                        out.text (".f.");
    
    out.text ("                       ");
    out.line (max_level_of_large_forces_);
}

void M_DynaMix_input::
write_RDFs (Text_output & out)
{
    out.line ("#");
    out.line ("#  RDFs have a separate restart file with extension .rdf");
    out.line ("#  Calculate RDF?      Dump RDF?     Read restart RDF file?");

    out.text ("  ");
    if (calculate_RDF_)         out.text (".t.");
    else                        out.text (".f.");

    out.text ("                  ");
    if (dump_RDF_)              out.text (".t.");
    else                        out.text (".f.");

    out.text ("           ");
    if (read_restart_RDF_file_) out.line (".t.");
    else                        out.line (".f.");

    out.line ("#");
    out.line ("#  RDF for all sites?");
    out.line ("#  (if .f., they should");
    out.line ("#  be specified below)    Cutoff-RDF(A)       Resolution of RDFs");

    out.text ("  ");
    if (is_RDF_for_all_sites_)  out.text (".t.");
    else                        out.text (".f.");

    out.text ("                    ");
    out.text (cutoff_RDF_);

    out.text ("                  ");
    out.line (RDF_resolution_);
}

void M_DynaMix_input::
write_dump_trajectory (Text_output & out)
{
    out.line ("#");
    out.line ("#  dump trajectory (0/1/2)?                number of config. in a trajectory");
    out.line ("#  0 - no trajectory                      file. Trajectory files have extensions");
    out.line ("#  1 - unformatted files                   .001, .002, ...");
    out.line ("#  2 - \".xmol\" format         interval(s)");

    out.text ("  ");
    out.text (dump_trajectory_format_);

    out.text ("                           ");
    out.text (dump_trajectory_interval_);

    out.text ("           ");
    out.line (number_of_config_in_a_trajectory_);
    
    out.line ("#  dump trajectory of molecules (1/0):");

    for (int i=0;  i<molecular_types_.size();  ++i)
    {
        out.text ("  ");
        out.text (molecular_types_[i].dump_trajectory_);
    }
    out.line ("");
}

void M_DynaMix_input::
write_TCF_calculations (Text_output & out)
{
    out.line ("#");
    out.line ("#   TCF calculations");
    out.line ("#   ----------------");
    out.line ("#  Attention!  Tcf calculations are NOT parallelized and may essentially");
    out.line ("#  slow down parallel simulations");
    out.line ("#");
    out.line ("#  TCFs have a separate restart file with extension .tcf");
    out.line ("#  NSTEG is number of points for calculation TCF");
    out.line ("#  JUMP - number of MD steps between the points for calculation of tcf.");
    out.line ("# calculate tcf?   restart tcf?   dump tsf?    NSTEG      JUMP");

    out.text ("  ");
    if (calculate_tcf_)     out.text (".t.");
    else                    out.text (".f.");

    out.text ("            ");
    if (restart_tcf_)       out.text (".t.");
    else                    out.text (".f.");

    out.text ("           ");
    if (dump_tsf_)          out.text (".t.");
    else                    out.text (".f.");
    
    out.text ("          ");
    out.text (NSTEG_);

    out.text ("        ");
    out.line (JUMP_);

    out.line ("#");
    out.line ("#  which of 12 tcf calculate (0/1/2):");
    out.line ("#  12 types:");
    out.line ("#  1 - velocity autocorellations");
    out.line ("#  2 - angular velocity autocorellations");
    out.line ("#  3 - 1 order Legendre polynom for dipole moment");
    out.line ("#  4 - 2 order Legendre polynom for dipole moment"); 
    out.line ("#  5 - 1 order Legendre polynom for reorintational tcf defined by vector below");
    out.line ("#  6 - 2 order Legendre polynom for reorintational tcf defined by vector below"); 
    out.line ("#  7 - X projection of velocity TCF");
    out.line ("#  8 - Y projection of velocity TCF");
    out.line ("#  9 - Z projection of velocity TCF");   
    out.line ("#  10 - X projection of angular velocity TCF");
    out.line ("#  11 - Y projection of angular velocity TCF");
    out.line ("#  12 - Z projection of angular velocity TCF");
    out.line ("#  if \"2\" specified for tcf 7-9 or 10-12, then tcf projections are");
    out.line ("#  calculated in molecular principal coordinate system; otherwize");
    out.line ("#  they calculated in laboratorial coordinate system");

    //fix
    out.line ("     1 1 1 1 1 1 2 2 2 0 0 0");
    
    out.line ("#");
    out.line ("#   unit vectors for reorientational tcf");
    out.line ("#   this vector is defined by 2 selected atoms on each molecule");

    //fix
    out.line ("    1    1    1    1    1    1    1    1    1    1");
    out.line ("    2    1    1    1    1    1    1    1    1    1");  
}

void M_DynaMix_input::
write_optional_parameters (Text_output & out)
{
    out.line ("#");
    out.line ("#  Other optional parameters");
    out.line ("#  if line \"add <n>\" is omitted, no optional parameters");
    out.line ("#  Description of optional parameters is given in file Extra_param");
    out.line ("#");
    out.line ("add       4");
    out.line ("no aver");
    out.line ("visual");
    out.line ("gather");
    out.line ("T-sep");
}

//fix
void M_DynaMix_input::
write_RDFs_additional (Text_output & out)
{
    out.line ("#  Total number of RDF-s");
    out.line ("0");
    /*out.line ("#  RDFs are defined by \"global\" site number");
    out.line ("#  First molecule type H2O have 3 sites with number 1,2,3");
    out.line ("#  Second molecule Na has one site with number 4");
    out.line ("#  and third molecule Cl has one site with number 5");
    out.line ("#   O - O");
    out.line ("#  This means RDF between sites 1 and 1, i.e. O atoms on H2O molecule");
    out.line ("1  1");
    out.line ("#   O - H");
    out.line ("#  symbol & followed by a number means that these RDF will be averaged");
    out.line ("#  (in this case: two hydrogens atoms (atoms 2 and 3 in H2O molecule)");
    out.line ("#   are equivalent)");
    out.line ("&2");
    out.line ("1  2");
    out.line ("1  3");
    out.line ("#    O - Na");
    out.line ("#  atoms on the next molecule type (Na) have site number 4, on Cl - 5.");
    out.line ("1  4");
    out.line ("#    O - Cl");
    out.line ("1  5");
    out.line ("#    H - H");
    out.line ("&3");
    out.line ("2  2");
    out.line ("2  3");
    out.line ("3  3");      
    out.line ("#     H - Na");
    out.line ("&2");
    out.line ("2  4");
    out.line ("3  4");        
    out.line ("#     H - Cl");
    out.line ("&2");
    out.line ("2  5");
    out.line ("3  5");        
    out.line ("#     Na - Na");
    out.line ("4  4");
    out.line ("#     Na - Cl");
    out.line ("4  5");
    out.line ("#   Cl - Cl");
    out.line ("5  5");//*/
}

bool M_DynaMix_input::
next_line (std::ifstream & in, std::string & line)
{
    while (std::getline (in, line))
    {
        ++line_counter_;
        if (line[0] != '#' &&                                    //not a comment
            line.find_first_not_of (" \t") != std::string::npos) //not a blank line
            return true;
    }
    return false;
}

void M_DynaMix_input::
load (const MM::Text & file_name)
{
    using namespace std;
    source_file_name_ = file_name;

    ifstream   file (source_file_name_.c_str());
    string     line;
    line_counter_ = 0;

    if (file.fail())
    {
        to_user(). error (Text("Cannot open file " + source_file_name_));
        return;
    }

    try
    {
        //log () << line.c_str() << "\n";
        {
        next_line (file, line);
        istringstream   stream (line);
        stream >> output_control_;
        if (stream.fail() || output_control_ < 2 || output_control_ > 10)
            throw Text (line.c_str()) + 
                "\nOutput control parameter should be in range 2-10.";
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        stream >> main_file_name_;
        if (stream.fail())
            throw Text (line.c_str()) + 
                "\nCan not read base file name for output files.";
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        stream >> molecular_database_path_;
        if (stream.fail())
            throw Text (line.c_str()) + 
                "\nCan not read molecular database path.";
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        string  read_from_restart_file, 
                dump_restart_file, 
                check_only, 
                zero_counter_of_cpu_time;
        
        stream >> read_from_restart_file 
               >> dump_restart_file 
               >> check_only 
               >> zero_counter_of_cpu_time;
        
        if (stream.fail() || 
            !(read_from_restart_file == ".t."   || read_from_restart_file == ".f.") ||
            !(dump_restart_file == ".t."        || dump_restart_file == ".f.") ||
            !(check_only == ".t."               || check_only == ".f.") ||
            !(zero_counter_of_cpu_time == ".t." || zero_counter_of_cpu_time == ".f."))
            throw Text (line.c_str()) + 
                "\nCan not read restart file options.";
        read_from_restart_file_   = read_from_restart_file   == ".t." ? true : false;
        dump_restart_file_        = dump_restart_file        == ".t." ? true : false;
        check_only_               = check_only               == ".t." ? true : false;
        zero_counter_of_cpu_time_ = zero_counter_of_cpu_time == ".t." ? true : false;
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        string  constant_temperature, 
                constant_pressure, 
                anizotropic_NPT;
        
        stream >> constant_temperature 
               >> constant_pressure 
               >> anizotropic_NPT;
        
        if (stream.fail() || 
            !(constant_temperature == ".t." || constant_temperature == ".f.") ||
            !(constant_pressure == ".t."    || constant_pressure == ".f.") ||
            !(anizotropic_NPT == ".t."      || anizotropic_NPT == ".f."))
            throw Text (line.c_str()) + 
                "\nCan not read type of statistical ensemble.";
        constant_temperature_   = constant_temperature  == ".t." ? true : false;
        constant_pressure_      = constant_pressure     == ".t." ? true : false;
        anizotropic_NPT_        = anizotropic_NPT       == ".t." ? true : false;
        }
        {
        next_line (file, line); // Number of molecular types 
        next_line (file, line); // Types of molecules 
        next_line (file, line); // Number of molecules 
        next_line (file, line); // non-bonded intramolecular interactions 
        next_line (file, line); // 1 - 4 intramolecular interactions
        next_line (file, line); // Scaling factors for 1-4 LJ interactions
        next_line (file, line); // Scaling factor for 1-4 electrostatic interactions
        next_line (file, line); // intramolecular potential type
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        stream >> temperature_ >> density_ >> pressure_;
        if (stream.fail())
            throw Text (line.c_str()) + "\nCan not read: "
                "temperature(K), density(g/cm^3), pressure(atm).";
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        double x, y, z;
        stream >> x >> y >> z >> cell_type_;
        if (stream.fail())
            throw Text (line.c_str()) + "\nCan not read: "
                "Rules for initial box size.";
        model().kit().boundary_conditions().box().set (x,y,z, 90,90,90);
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        //stream >> time_step_ >> small_steps_in_one_long_;
        //time_step_ *= 1e15;
        double time_step;
        stream >> time_step >> small_steps_in_one_long_;
        time_step_.set_fs (time_step);

        if (stream.fail())
            throw Text (line.c_str()) + "\nCan not read: "
                "Time step (s), Small steps in one long.";
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        //stream >> total_long_steps_ 
        int total_long_steps;
        stream >> total_long_steps 
               >> steps_for_intermediate_averaging_ 
               >> take_averages_each___steps_
               >> dump_restart_file_after___steps_;

        if (stream.fail())
            throw Text (line.c_str()) + "\nCan not read: "
                "Total (long) MD-steps, "
                "Steps for intermediate averaging, "
                "take averages each .. steps, "
                "dump restart file after .. steps.";

        time_step_.set_steps (total_long_steps);
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        string is_simple_velosity_scaling;
        stream >> t_termostat_ 
               >> p_termostat_ 
               >> is_simple_velosity_scaling
               >> delta_t_;
        if (stream.fail() ||
            !(is_simple_velosity_scaling == ".t." || 
              is_simple_velosity_scaling == ".f."))
            throw Text (line.c_str()) + "\nCan not read: "
                "T-termostat param(fs), "
                "P-termostat(fs), "
                "Simple velosity scaling, "
                "delta T (K)";
        is_simple_velosity_scaling_ = is_simple_velosity_scaling == ".t." 
            ? true : false;
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        stream >> r_cutoff_ 
               >> r_cutoff_fast_force_ 
               >> check_neigbours_after___step_;
        if (stream.fail())
            throw Text (line.c_str()) + "\nCan not read: "
                "Rcutoff(A), Rcut-fast forces, check neigbours after .. steps";
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        stream >> alpha_r_ 
               >> fexp_;
        if (stream.fail())
            throw Text (line.c_str()) + "\nCan not read: "
                "alpha*R, fexp (m/s**2)";
        }
        {
        next_line (file, line); // which of the molecules move 
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        string  recalculate_list_of_intramolecular_interactions;
        
        stream >> recalculate_list_of_intramolecular_interactions; 
        
        if (stream.fail() || 
            !(recalculate_list_of_intramolecular_interactions == ".t." || 
              recalculate_list_of_intramolecular_interactions == ".f."))
            throw Text (line.c_str()) + "\nCan not read: "
                "recalculate the list of intramolecular interactions.";
        recalculate_list_of_intramolecular_interactions_ = 
            recalculate_list_of_intramolecular_interactions  == ".t." 
                ? true : false;
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        string  constrain_dynamics;
        
        stream >> constrain_dynamics >> tolerance_parameter_; 
        
        if (stream.fail() || 
            !(constrain_dynamics == ".t." || constrain_dynamics == ".f."))
            throw Text (line.c_str()) + "\nCan not read: "
                "Constrain dynamics, tolerance parameter.";
        constrain_dynamics_ = constrain_dynamics  == ".t." ? true : false;
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        string  initial_state, set_velocities_to_0;
        
        stream >> initial_state >> set_velocities_to_0; 
        
        if (stream.fail() || 
            !(set_velocities_to_0 == ".t." || set_velocities_to_0 == ".f."))
            throw Text (line.c_str()) + "\nCan not read: "
                "set velocities to 0.";
        set_velocities_to_0_ = set_velocities_to_0  == ".t." ? true : false;
        }
        {
        next_line (file, line); // Parameters for above in cases 2,3  
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        string  keep_initial_configuration, file_name;
        
        stream >> hole_radius_ 
               >> keep_initial_configuration 
               >> permissible_deviation_ 
               >> file_name; 
        
        if (stream.fail() || 
            !(keep_initial_configuration == ".t." || 
              keep_initial_configuration == ".f."))
            throw Text (line.c_str()) + "\nCan not read: "
                "Radius of a hole, Keep initial configuration, "
                "Permissible deviation, File name.";
        keep_initial_configuration_ = keep_initial_configuration  == ".t." 
            ? true : false;
        file_name_ = file_name.c_str();
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        string  change_temperature_at_restart, change_density_at_restart;
        
        stream >> change_temperature_at_restart 
               >> change_density_at_restart; 
        
        if (stream.fail() || 
            !(change_temperature_at_restart == ".t." || 
              change_temperature_at_restart == ".f." || 
              change_density_at_restart     == ".t." || 
              change_density_at_restart     == ".f."))
            throw Text (line.c_str()) + "\nCan not read: "
                "Change temperature at restart,  change density.";
        change_temperature_at_restart_ = change_temperature_at_restart == ".t."
            ? true : false;
        change_density_at_restart_ = change_density_at_restart == ".t."
            ? true : false;
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        string  dump_XMOL_config_file;
        
        stream >> final_averaging_after___intermediate_averaging_ 
               >> dump_XMOL_config_file; 
        
        if (stream.fail() || 
            !(dump_XMOL_config_file == ".t." || 
              dump_XMOL_config_file == ".f."))
            throw Text (line.c_str()) + "\nCan not read: "
                "Final averaging after .. intermediate averaging, "
                "Dump XMOL config. file.";
        dump_XMOL_config_file_ = dump_XMOL_config_file == ".t." ? true : false;
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        string  cut_large_forces;
        
        stream >> cut_large_forces 
               >> max_level_of_large_forces_; 
        
        if (stream.fail() || 
            !(cut_large_forces == ".t." || 
              cut_large_forces == ".f."))
            throw Text (line.c_str()) + "\nCan not read: "
                "Cut large forces, Parameter";
        cut_large_forces_ = cut_large_forces == ".t." ? true : false;
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        string  calculate_RDF, dump_RDF, read_restart_RDF_file;
        
        stream >> calculate_RDF 
               >> read_restart_RDF_file 
               >> dump_RDF; 
        
        if (stream.fail() || 
            !(calculate_RDF == ".t." || 
              calculate_RDF == ".f." ||
              dump_RDF      == ".t." || 
              dump_RDF      == ".f." ||
              read_restart_RDF_file      == ".t." || 
              read_restart_RDF_file      == ".f."))
            throw Text (line.c_str()) + "\nCan not read: "
                "Calculate RDF?  Dump RDF? Read restart RDF file";
        calculate_RDF_         = calculate_RDF         == ".t." ? true : false;
        dump_RDF_              = dump_RDF              == ".t." ? true : false;
        read_restart_RDF_file_ = read_restart_RDF_file == ".t." ? true : false;
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        string  is_RDF_for_all_sites;
        
        stream >> is_RDF_for_all_sites 
               >> cutoff_RDF_ 
               >> RDF_resolution_; 
        
        if (stream.fail() || 
            !(is_RDF_for_all_sites == ".t." || 
              is_RDF_for_all_sites == ".f."))
            throw Text (line.c_str()) + "\nCan not read: "
                "RDF for all sites?  Cutoff-RDF(A),  Resolution of RDFs.";
        is_RDF_for_all_sites_  = is_RDF_for_all_sites == ".t." ? true : false;
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        
        stream >> dump_trajectory_format_ 
               >> dump_trajectory_interval_ 
               >> number_of_config_in_a_trajectory_; 
        
        if (stream.fail()) 
            throw Text (line.c_str()) + "\nCan not read: "
                "dump trajectory,  interval(s),  number of config.";
        if (!(dump_trajectory_format_ == 0 ||
              dump_trajectory_format_ == 1 ||
              dump_trajectory_format_ == 2 ))
            throw Text (line.c_str())+"\nDump trajectory should be 0, 1 or 2.";
        }
        {
        next_line (file, line); // dump trajectory of molecules (1/0):  
        }
        {
        next_line (file, line);
        istringstream   stream (line);
        string  calculate_tcf, restart_tcf, dump_tsf;
        
        stream >> calculate_tcf 
               >> restart_tcf 
               >> dump_tsf 
               >> NSTEG_ 
               >> JUMP_; 
        
        if (stream.fail() || 
            !(calculate_tcf == ".t." || 
              calculate_tcf == ".f." ||
              restart_tcf   == ".t." || 
              restart_tcf   == ".f." ||
              dump_tsf      == ".t." || 
              dump_tsf      == ".f."))
            throw Text (line.c_str()) + "\nCan not read: "
                "calculate tcf?   restart tcf?   dump tsf?,  NSTEG,  JUMP";
        calculate_tcf_ = calculate_tcf == ".t." ? true : false;
        restart_tcf_   = restart_tcf   == ".t." ? true : false;
        dump_tsf_      = dump_tsf      == ".t." ? true : false;
        }
    }
    catch (Text & error)
    {
        to_user().error (Text() + 
            "File " + source_file_name_ + "\n"
            "is broken at line " + line_counter_ + ":\n" +
            error);
    }
}

}//MM







