diff --git a/Documentation/source/guides/forcefield.rst b/Documentation/source/guides/forcefield.rst index bc9b3066..dd2feb82 100644 --- a/Documentation/source/guides/forcefield.rst +++ b/Documentation/source/guides/forcefield.rst @@ -244,6 +244,37 @@ combining rules can be used, as described elsewhere. Energy parameter :math:`\epsilon/k_B` K Collision diameter :math:`\sigma` Å =================== ====================== ================ + +Mixing rules +~~~~~~~~~~~~ + +The available mixing rules are Lorentz-Berthelot, geometric, size-weighted and custom. +The algebraic expressions for the first three are: + +- Lorenz-Berthelot + + .. math:: + \epsilon_{ij} = \sqrt {\epsilon_i \epsilon_j} \\ + \sigma_{ij} = \frac{\sigma_i + \sigma_j}{2} + +- Geometric + + .. math:: + \epsilon_{ij} = \sqrt {\epsilon_i \epsilon_j} \\ + \sigma_{ij} = \sqrt{\sigma_i \sigma_j} + +- Size-weighted + + .. math:: + \epsilon_{ij} = \frac{ \frac{\epsilon_i}{\epsilon_j} \sigma_i + \frac{\epsilon_j}{\epsilon_i} \sigma_j }{2} \\ + \sigma_{ij} = \frac{\sigma_i + \sigma_j}{2} + +where :math:`\epsilon_{ij}` and :math:`\sigma_{ij}` the energy and size parameters +for atom types `i` and `j`, respectively. + +In addition to these analytic mixing rules, the user can specify the specific +interaction parameters for each pair of atom types using the ``custom`` +mixing rule. For more information, see :ref:`sec:mixing_rule` Electrostatics ~~~~~~~~~~~~~~ diff --git a/Documentation/source/guides/input_files.rst b/Documentation/source/guides/input_files.rst index a33a4d2e..ff65bcfc 100644 --- a/Documentation/source/guides/input_files.rst +++ b/Documentation/source/guides/input_files.rst @@ -44,6 +44,8 @@ For example:: | Cassandra will then use ``dee.out`` as prefix for all output files created. +.. _sec:sim_type: + Simulation Type ~~~~~~~~~~~~~~~ @@ -79,6 +81,8 @@ temperature :math:`T` are held constant. Simulation types ``fragment`` and ``ring_fragment`` are used only for generating a fragment library. +.. _sec:num_species: + Number of species ~~~~~~~~~~~~~~~~~ @@ -197,6 +201,8 @@ selection of *Character(i,2)* as described below: For all options, cutoff distances must be less than or equal to the shortest edge length of a simulation box. +.. _sec:charge_style: + Charge Style ~~~~~~~~~~~~ @@ -291,6 +297,8 @@ liquid phase simulations. the half of the simulation box length. In the case of Ewald summation, however, the accuracy will be the same as *Real(i,4)*. +.. _sec:mixing_rule: + Mixing Rule ~~~~~~~~~~~ @@ -299,9 +307,13 @@ Mixing Rule Sets the method by which van der Waals interactions between unlike atoms are calculated. Acceptable options are ``lb`` for Lorentz-Berthelot, ``geometric`` -for geometric mixing rule and ``custom`` for allowing the user to provide -specific values. To use either ``lb`` or ``geometric`` keywords with the Mie +for geometric mixing rule, ``sw`` for size-weighted +and ``custom`` for allowing +the user to provide specific values. To use either ``lb`` or ``geometric`` +keywords with the Mie potential, all atomtypes must have the same repulsive and dispersive exponents. +Currently, +the size-weighted mixing rule is not implemented for the Mie potential. If this section is missing, ``lb`` is used as default. To illustrate the use of the ``custom`` option, consider a mixture of methane @@ -360,6 +372,7 @@ As an example, checkpoint file. If ``# Start_Type`` is ``read_config``, then the seeds specified in the input file are used. +.. _sec:min_cutoff: Minimum Cutoff ~~~~~~~~~~~~~~ @@ -379,6 +392,8 @@ important that the minimum distance is set to be less than the shortest distance between any two sites on the molecule. For most systems, 1 Å seems to work OK, but for models with dummy sites, a shorter value may be required. +.. _sec:pair_storage: + Pair Energy Storage ~~~~~~~~~~~~~~~~~~~ @@ -391,6 +406,8 @@ thereby saving a loop over all molecules. This requires more memory, but it can be faster. The default is to not use this feature. If you wish to use this, set *Logical* to ``true``. +.. _sec:mol_files: + Molecule Files ~~~~~~~~~~~~~~ @@ -421,6 +438,8 @@ molecules and 5 octane molecules in the total system. The maximum number of molecules specified here will be used to allocate memory for each species, so do not use larger numbers than are needed. +.. _sec:sim_box: + Simulation Box ~~~~~~~~~~~~~~ @@ -496,12 +515,12 @@ several options: ``sphere``, ``cylinder``, ``slitpore``, or ``interface``. Each option requires additional parameters, as follows: - | ``sphere r``, where ``r`` is the radius of a sphere centered at the origin + - | ``cylinder r``, where ``r`` is the radius of a cylinder centered on the z-axis -- | ``slitpore z_max``, where ``z_max`` is half the height of a rectangular prism - centered on the *xy*-plane -- | ``interface z_min z_max``, which defines two rectangular prisms that span - the box in the *x* and *y* directions. One box has bounds ``z_min < z < z_max`` - and the other has bounds ``-z_max < z < -z_min``. + +- | ``slitpore z_max``, where ``z_max`` is half the height of a rectangular prism centered on the *xy*-plane + +- | ``interface z_min z_max``, which defines two rectangular prisms that span the box in the *x* and *y* directions. One box has bounds ``z_min < z < z_max`` and the other has bounds ``-z_max < z < -z_min``. For example, to make a spherical droplet with a radius of 5 Å in cubic box with 100 Å side lengths: @@ -517,6 +536,7 @@ with 100 Å side lengths: In addition, the insertion method for each species must be identified in the ``Start_Type`` or ``Move_Probability_Info`` sections. +.. _sec:temperature: Temperature ~~~~~~~~~~~ @@ -533,6 +553,8 @@ the temperature of box 2 will be read from a second line: 300.0 300.0 +.. _sec:pressure: + Pressure ~~~~~~~~ @@ -551,6 +573,8 @@ the pressure of box 2 will be read from a second line: If the simulation type does not require an input pressure (e.g., NVT), this section will be ignored. +.. _sec:chemical_potential: + Chemical Potential ~~~~~~~~~~~~~~~~~~ @@ -584,6 +608,8 @@ potential (for methane) is required and the following are equivalent: Specifying the chemical potential as ``0.0`` is **not** the same as ``none``. +.. _sec:move_probabilities: + Move Probabilities ~~~~~~~~~~~~~~~~~~ @@ -602,6 +628,8 @@ subsections must be between the section header and footer. If the move probabilities do not sum to 1.0, then the probability of each move will be divided by the total. +.. _sec:translation: + Translation ^^^^^^^^^^^ @@ -633,6 +661,8 @@ displacement of the solid species to zero. Every molecule in the simulation with a maximum displacement greater than zero has an equal chance of being moved. +.. _sec:rotation: + Rotation ^^^^^^^^ @@ -686,6 +716,8 @@ Molecules of species 1 will be rotated a maximum of 30 around the x, y or z-axis, molecules of species 2 will be rotated by choosing a random unit vector, and the point particles will not be rotated. +.. _sec:angle: + Angle ^^^^^ @@ -714,6 +746,7 @@ containing angles within a given box. the case where parts of its fragments are rarely regrown by a regrowth move. +.. _sec:dihedral: Dihedral ^^^^^^^^ @@ -751,6 +784,8 @@ systems containing more than one box. case where the parts of its fragments are rarely regrown (albeit a small maximum width is provided). +.. _sec:regrowth: + Regrowth ^^^^^^^^ @@ -786,6 +821,8 @@ balance of regrowth moves on a molecule of species 2. *Real(2,i)* should be set to zero for monatomic, linear, or rigid species, including solid frameworks. +.. _sec:volume: + Volume ^^^^^^ @@ -841,6 +878,8 @@ achieve about 50% acceptance rates. The volume perturbation move is only supported for cubic boxes. +.. _sec:insert_delete: + Insertion and Deletion Moves ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -961,6 +1000,8 @@ donor box for 75% of attempted swaps, use: The probability of selecting a species with insertion method ``none`` must be 0. +.. _sec:ring_flip: + Ring Flip Move ^^^^^^^^^^^^^^ | ``# Prob_Ring`` @@ -1181,6 +1222,7 @@ rejected and re-attempted. read_config 100 0 water.xyz make_config  50  50 +.. _sec:run_type: Run Type ~~~~~~~~ @@ -1417,6 +1459,7 @@ This results in a dramatic speed boost for Widom insertions in dense media such most trial positions have core overlap. This method only supports cubic boxes and orthogonal boxes. +.. _sec:prop_output: Property Output ~~~~~~~~~~~~~~~ @@ -1503,6 +1546,8 @@ may be output using the following format: density pressure +.. _sec:fragment_files: + Fragment Files ~~~~~~~~~~~~~~ @@ -1537,6 +1582,8 @@ This tells Cassandra to use the files ``frag_1_1.dat``, ``frag_2_1.dat`` and identical fragment, Cassandra will use the same fragment library ``frag_1_2.dat`` for these fragments. +.. _sec:verbosity_logfile: + Verbosity in log file ~~~~~~~~~~~~~~~~~~~~~ @@ -1549,6 +1596,8 @@ for development purposes. If this section is missing, *Logical* is set to ``false`` by default. Supported options for *Logical* are ``true`` or ``false``. +.. _sec:file_info: + File Info ~~~~~~~~~ @@ -1572,6 +1621,8 @@ the fragment being simulated. This tells Cassandra to store the fragment library in the file named ``frag.dat``. +.. _sec:cbmc_parameters: + CBMC parameters ~~~~~~~~~~~~~~~ diff --git a/Src/create_nonbond_table.f90 b/Src/create_nonbond_table.f90 index 42457ec4..3dd9c9f8 100755 --- a/Src/create_nonbond_table.f90 +++ b/Src/create_nonbond_table.f90 @@ -426,6 +426,64 @@ SUBROUTINE Create_Nonbond_Table ENDDO ENDDO + + ELSE IF (mix_rule == 'SW') THEN + + DO itype = 1, nbr_atomtypes + DO jtype = itype+1, nbr_atomtypes + + IF (int_vdw_style(1) == vdw_lj) THEN + ! There are two vdw parameters that need mixing + + ! epsilonij needs both epsilons and sigmas + + eps_i = vdw_param1_table(itype,itype) + eps_j = vdw_param1_table(jtype,jtype) + + sig_i = vdw_param2_table(itype,itype) + sig_j = vdw_param2_table(jtype,jtype) + + ! for parameters with zero, avoid overflow and set to zero + IF (eps_i <= tiny_number) eps_i = 0.0_DP + IF (eps_j <= tiny_number) eps_j = 0.0_DP + + IF (sig_i <= tiny_number .OR. sig_j <= tiny_number) THEN + err_msg = "" + err_msg(1) = "The size-weighted mixing rule requires a non-zero sigma LJ parameter" + CALL Clean_Abort(err_msg,'Create_Nonbond_Table') + END IF + + ! SW mixing_rule: epsij = (sigma_i / sigma_j) * epsilon_i + (sigma_j / sigma_i) * epsilon_j + vdw_param1_table(itype,jtype) = ( (sig_i / sig_j) * eps_i + (sig_j / sig_i) * eps_j ) * 0.5_DP + + ! SW mixing rule: sigmaij = (sigmai * sigmaj)^(1/2) + vdw_param2_table(itype,jtype) = (sig_i + sig_j) * 0.5_DP + + ! Report parameters to logfile. + IF (verbose_log) THEN + WRITE(logunit,'(X,A23,5X,A23,2X,F12.4,X,F12.4)') & + atom_type_list(itype), atom_type_list(jtype), & + vdw_param1_table(itype,jtype), vdw_param2_table(itype,jtype) + ENDIF + + ! Mixed interactions are symmetric + vdw_param1_table(jtype,itype) = vdw_param1_table(itype,jtype) + vdw_param2_table(jtype,itype) = vdw_param2_table(itype,jtype) + + ! Mark that the parameters have been set + vdw_param_set(itype,jtype) = 1 + vdw_param_set(jtype,itype) = 1 + + ELSEIF (int_vdw_style(1) == vdw_mie) THEN + err_msg = "" + err_msg(1) = "The size-weighted mixing rule has not been implemented for the Mie potential yet." + CALL Clean_Abort(err_msg,'Create_Nonbond_Table') + END IF + + ENDDO + ENDDO + + ELSE IF (mix_rule == 'custom') THEN REWIND(inputunit) diff --git a/Src/input_routines.f90 b/Src/input_routines.f90 index 835aba90..24e4b1ba 100755 --- a/Src/input_routines.f90 +++ b/Src/input_routines.f90 @@ -936,6 +936,9 @@ SUBROUTINE Get_Mixing_Rules WRITE(logunit,'(A)') 'Lorentz-Berthelot mixing rule specified' ELSEIF (mix_rule == 'geometric') THEN WRITE(logunit,'(A)') 'Geometric mixing rule specified' + ELSEIF (mix_rule == 'sw'.OR. mix_rule == 'SW') THEN + mix_rule = 'SW' + WRITE(logunit,'(A)') 'Size-weighted mixing rule specified' ELSEIF (mix_rule == 'custom') THEN WRITE(logunit,'(A)') 'Custom mixing rule specified' ELSE