diff --git a/GUIDE.md b/GUIDE.md index 500922f..c6ac74e 100644 --- a/GUIDE.md +++ b/GUIDE.md @@ -789,7 +789,7 @@ composed of Lennard-Jones interaction sites. For simplicity the sites are taken to be identical, although the program is easily generalized. Molecular orientations are represented by quaternions, which are used to calculate the rotation matrix -and hence the LJ site positions. +and hence the interaction site positions. We test this with the three-site model of orthoterphenyl, a fragile glassformer, described in the following publications amongst others. @@ -799,46 +799,54 @@ described in the following publications amongst others. * S Mossa, E La Nave, HE Stanley, C Donati, F Sciortino, P Tartaglia, _Phys Rev E,_ __65,__ 041205 (2002) * E La Nave, S Mossa, F Sciortino, P Tartaglia, _J Chem Phys,_ __120,__ 6128 (2004) +We compare with the results of Mossa et al (2002). The sites are arranged at the vertices of an isosceles triangle with bond angle 75 degrees, -LJ parameters ε/_k_B = 600K -or ε = 5 kJ mol-1 to a good approximation, +LJ parameters ε = 5.276 kJ mol-1, σ=0.483nm, -and bond lengths equal to σ. +and two equal bonds of length σ. +The program employs the usual reduced units based on ε and σ +and in these units the potential cutoff of Mossa et al (2002) is _R_c=2.612; +the pair potential is Lennard-Jones with a shifted-force correction term, linear in _r_, +to make the potential and its derivative vanish at _r_=_R_c. +Apart from the temperatures, +default program parameters were used throughout the tests. Tests were performed at ρ=0.32655 which is equivalent to ρ4=1.108g cm-3 in Mossa et al (2002). -Comparisons of potential energy (_PE_=_E_-3 _T_ converted to kJ/mol with a factor 5) -were made with the fit given by eqn (6) of that paper. - -_T_ | _E_ | _P_ | _T_ (K) | _PE_ (kJ/mol) | eqn (6) ------ | ----- | ----- | ----- | ----- | ----- -0.5 | -14.30(1) | 1.63(3) | 300 | -79.00(5) | -77.52 -1.0 | -11.21(1) | 5.52(3) | 600 | -71.05(5) | -68.55 -1.5 | -8.262(7) | 8.95(2) | 900 | -63.81(4) | -61.19 -2.0 | -5.52(1) | 11.86(3) | 1200 | -57.60(5) | -54.69 - -Exact agreement is not expected because the potential of Mossa et al (2002) has a different -cutoff correction, but the agreement is reasonable. -A second set of tests were performed at _T_=0.63333=380K -at the specified densities ρ1, … ρ5. +Comparisons of potential energy (_PE_=_E_-3 _T_ converted to kJ/mol with a factor 5.276) +were made with the fit given by eqn (23) of that paper. +Note that ε/kB≈635 K. + +_T_ | _E_ | _P_ | _T_ (K) | _PE_ (kJ/mol) | _PE_ (kJ/mol) eqn (23) +----- | ----- | ----- | ----- | ----- | ----- +0.5 | -12.993(3) | 1.773(8) | 317 | -76.47(2) | -76.945 +1.0 | -9.817(9) | 5.86(2) | 635 | -67.62(5) | -67.634 +1.5 | -6.84(1) | 9.38(4) | 952 | -59.83(5) | -60.011 +2.0 | -4.07(1) | 12.37(4) | 1270 | -53.13(5) | -53.265 + +A second set of tests was performed at _T_=0.6≈380K +at the specified densities ρ1, … ρ5 of Mossa et al (2002). Here the excess pressure (_P_(ex)=_P_-ρ_T_ converted to MPa -with a factor 73.54 based on the values of ε and σ) +with a factor 77.75 based on the values of ε and σ) is compared with the fit given by eqn (28) and the coefficients in Table III of Mossa et al (2002). NB the volumes to insert into the equation are those of their Table I, which are specific to their system size. - -Id | ρ | _E_ | _P_ | _P_(ex) (MPa) | eqn (28) ------ | ----- | ----- | ----- | ----- | ----- -1 | 0.30533 | -12.737(7) | 0.35(2) | 12(2) | 19.077 -2 | 0.31240 | -13.025(7) | 0.99(2) | 58(2) | 60.143 -3 | 0.31918 | -13.293(8) | 1.66(2) | 107(2) | 112.798 -4 | 0.32655 | -13.50(1) | 2.60(2) | 176(1) | 177.222 -5 | 0.33451 | -13.65(1) | 3.95(2) | 274(1) | 253.510 - -Although not perfect at the ends of the range, the agreement is not bad; -once more, the difference in cutoff correction should be borne in mind. -Note that Mossa et al (2002) used a different value -ε=5.276 kJ/mol which, with their cutoff term, gives a well depth 4.985 kJ/mol. +In addition their eqn (29) with coefficients in Table V is a fit to their potential energy, +which we calculate from the simulation as described above. + +Id | ρ | _E_ | _P_ | _P_(ex) (MPa) | _P_(ex) (MPa) eqn (28) | _PE_ (kJ/mol) | _PE_ (kJ/mol) eqn (29) +----- | ----- | ----- | ----- | ----- | ----- | ------ | ------ +1 | 0.30533 | -11.625(4) | 0.45(1) | 20.7(7) | 19.077 | -70.83(2) | -70.818 +2 | 0.31240 | -11.914(6) | 1.01(2) | 64(2) | 60.143 | -72.36(3) | -72.289 +3 | 0.31918 | -12.143(5) | 1.74(1) | 120(2) | 112.798 | -73.56(3) | -73.601 +4 | 0.32655 | -12.400(4) | 2.53(1) | 181(1) | 177.222 | -74.92(2) | -74.886 +5 | 0.33451 | -12.487(4) | 3.84(1) | 283(1) | 253.510 | -75.38(2) | -75.825 + +In making these comparisons, +our default run length (10 blocks of 20000 sweeps each) should be borne in mind, +since this system can show sluggish behaviour. +The MD simulations of Mossa et al (2002) are reported to extend to several hundred nanoseconds +(of order 107 MD timesteps) at the lowest temperatures. ## DPD program For the `dpd` example, we recommend generating an initial configuration diff --git a/mc_nvt_poly_lj.f90 b/mc_nvt_poly_lj.f90 index d54f079..b16d601 100644 --- a/mc_nvt_poly_lj.f90 +++ b/mc_nvt_poly_lj.f90 @@ -39,14 +39,14 @@ PROGRAM mc_nvt_poly_lj ! For example, for Lennard-Jones, sigma = 1, epsilon = 1 ! Despite the program name, there is nothing here specific to Lennard-Jones - ! The model is defined in mc_module + ! The model (including the cutoff distance) is defined in mc_module USE, INTRINSIC :: iso_fortran_env, ONLY : input_unit, output_unit, error_unit, iostat_end, iostat_eor USE config_io_module, ONLY : read_cnf_mols, write_cnf_mols USE averages_module, ONLY : run_begin, run_end, blk_begin, blk_end, blk_add - USE maths_module, ONLY : metropolis, random_rotate_quaternion, random_translate_vector + USE maths_module, ONLY : metropolis, random_rotate_quaternion, random_translate_vector, q_to_a USE mc_module, ONLY : introduction, conclusion, allocate_arrays, deallocate_arrays, & - & potential_1, potential, n, r, e, potential_type + & potential_1, potential, n, na, db, r, e, d, potential_type IMPLICIT NONE @@ -55,35 +55,34 @@ PROGRAM mc_nvt_poly_lj REAL :: dr_max ! Maximum MC translational displacement REAL :: de_max ! Maximum MC rotational displacement REAL :: temperature ! Specified temperature - REAL :: r_cut ! Potential cutoff distance ! Composite interaction = pot & vir & ovr variables TYPE(potential_type) :: total, partial_old, partial_new - INTEGER :: blk, stp, i, nstep, nblock, moves, ioerr - REAL :: delta, m_ratio - REAL, DIMENSION(3) :: ri - REAL, DIMENSION(0:3) :: ei + INTEGER :: blk, stp, i, a, nstep, nblock, moves, ioerr + REAL :: delta, m_ratio + REAL, DIMENSION(3) :: ri + REAL, DIMENSION(0:3) :: ei + REAL, DIMENSION(3,3) :: ai + REAL, DIMENSION(na,3) :: di CHARACTER(len=4), PARAMETER :: cnf_prefix = 'cnf.' CHARACTER(len=3), PARAMETER :: inp_tag = 'inp' CHARACTER(len=3), PARAMETER :: out_tag = 'out' CHARACTER(len=3) :: sav_tag = 'sav' ! May be overwritten with block number - NAMELIST /nml/ nblock, nstep, temperature, r_cut, dr_max, de_max + NAMELIST /nml/ nblock, nstep, temperature, dr_max, de_max WRITE ( unit=output_unit, fmt='(a)' ) 'mc_nvt_poly_lj' WRITE ( unit=output_unit, fmt='(a)' ) 'Monte Carlo, constant-NVT ensemble, polyatomic molecule' - WRITE ( unit=output_unit, fmt='(a)' ) 'Simulation uses cut-and-shifted potential' CALL introduction CALL RANDOM_SEED () ! Initialize random number generator ! Set sensible default run parameters for testing nblock = 10 - nstep = 10000 + nstep = 20000 temperature = 1.0 - r_cut = 2.5 dr_max = 0.05 de_max = 0.05 @@ -101,7 +100,6 @@ PROGRAM mc_nvt_poly_lj WRITE ( unit=output_unit, fmt='(a,t40,i15)' ) 'Number of blocks', nblock WRITE ( unit=output_unit, fmt='(a,t40,i15)' ) 'Number of steps per block', nstep WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Temperature', temperature - WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Potential cutoff distance', r_cut WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Maximum r displacement', dr_max WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Maximum e displacement', de_max @@ -110,13 +108,21 @@ PROGRAM mc_nvt_poly_lj WRITE ( unit=output_unit, fmt='(a,t40,i15)' ) 'Number of molecules', n WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Simulation box length', box WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Density of molecules', REAL(n) / box**3 - CALL allocate_arrays ( box, r_cut ) + CALL allocate_arrays ( box ) CALL read_cnf_mols ( cnf_prefix//inp_tag, n, box, r, e ) ! Second call is to get r and e r(:,:) = r(:,:) / box ! Convert positions to box units r(:,:) = r(:,:) - ANINT ( r(:,:) ) ! Periodic boundaries + ! Calculate all bond vectors + DO i = 1, n ! Loop over molecules + ai = q_to_a ( e(:,i) ) ! Rotation matrix for i + DO a = 1, na ! Loop over all atoms + d(:,a,i) = MATMUL ( db(:,a), ai ) ! NB: equivalent to ai_T*db, ai_T=transpose of ai + END DO ! End loop over all atoms + END DO ! End loop over molecules + ! Initial energy and overlap check - total = potential ( box, r_cut ) + total = potential ( box ) IF ( total%ovr ) THEN WRITE ( unit=error_unit, fmt='(a)') 'Overlap in initial configuration' STOP 'Error in mc_nvt_poly_lj' @@ -136,7 +142,7 @@ PROGRAM mc_nvt_poly_lj DO i = 1, n ! Begin loop over atoms - partial_old = potential_1 ( r(:,i), e(:,i), i, box, r_cut ) ! Old molecule potential, virial, etc + partial_old = potential_1 ( r(:,i), d(:,:,i), i, box ) ! Old molecule potential, virial, etc IF ( partial_old%ovr ) THEN ! should never happen WRITE ( unit=error_unit, fmt='(a)') 'Overlap in current configuration' @@ -146,8 +152,12 @@ PROGRAM mc_nvt_poly_lj ri = random_translate_vector ( dr_max/box, r(:,i) ) ! Trial move to new position (in box=1 units) ri = ri - ANINT ( ri ) ! Periodic boundary correction ei = random_rotate_quaternion ( de_max, e(:,i) ) ! Trial rotation + ai = q_to_a ( ei ) ! Rotation matrix for i + DO a = 1, na ! Loop over all atoms + di(:,a) = MATMUL ( db(:,a), ai ) ! NB: equivalent to ai_T*db, ai_T=transpose of ai + END DO ! End loop over all atoms - partial_new = potential_1 ( ri, ei, i, box, r_cut ) ! New molecule potential, virial etc + partial_new = potential_1 ( ri, di, i, box ) ! New molecule potential, virial etc IF ( .NOT. partial_new%ovr ) THEN ! Test for non-overlapping configuration @@ -155,10 +165,11 @@ PROGRAM mc_nvt_poly_lj delta = delta / temperature ! Divide by temperature IF ( metropolis ( delta ) ) THEN ! Accept Metropolis test - total = total + partial_new - partial_old ! Update potential energy - r(:,i) = ri ! Update position - e(:,i) = ei ! Update quaternion - moves = moves + 1 ! Increment move counter + total = total + partial_new - partial_old ! Update potential energy + r(:,i) = ri ! Update position + e(:,i) = ei ! Update quaternion + d(:,:,i) = di ! Update bond vectors + moves = moves + 1 ! Increment move counter END IF ! End accept Metropolis test END IF ! End test for non-overlapping configuration @@ -195,12 +206,11 @@ FUNCTION calc_variables ( ) RESULT ( variables ) ! This routine calculates all variables of interest and (optionally) writes them out ! They are collected together in the variables array, for use in the main program - ! In this example we simulate using the cut-and-shifted potential only - ! The values of < p_s >, < e_s > and density should be consistent (for this potential) + ! In this example we simulate using the shifted-force potential only + ! The values of < p_sf >, < e_sf > and density should be consistent (for this potential) ! There are no long-range or delta corrections - ! The value of the cut (but not shifted) potential is not used, in this example - TYPE(variable_type) :: m_r, e_s, p_s + TYPE(variable_type) :: m_r, e_sf, p_sf REAL :: vol, rho ! Preliminary calculations @@ -218,16 +228,16 @@ FUNCTION calc_variables ( ) RESULT ( variables ) ! Move acceptance ratio m_r = variable_type ( nam = 'Move ratio', val = m_ratio, instant = .FALSE. ) - ! Internal energy per molecule (cut-and-shifted potential) - ! Ideal gas contribution (assuming nonlinear molecules) plus total (cut-and-shifted) PE divided by N - e_s = variable_type ( nam = 'E/N cut&shift', val = 3.0*temperature + total%pot/REAL(n) ) + ! Internal energy per molecule (shifted-force potential) + ! Ideal gas contribution (assuming nonlinear molecules) plus total PE divided by N + e_sf = variable_type ( nam = 'E/N shifted force', val = 3.0*temperature + total%pot/REAL(n) ) - ! Pressure (cut-and-shifted potential) + ! Pressure (shifted-force potential) ! Ideal gas contribution plus total virial divided by V - p_s = variable_type ( nam = 'P cut&shift', val = rho*temperature + total%vir/vol ) + p_sf = variable_type ( nam = 'P shifted force', val = rho*temperature + total%vir/vol ) ! Collect together for averaging - variables = [ m_r, e_s, p_s ] + variables = [ m_r, e_sf, p_sf ] END FUNCTION calc_variables diff --git a/mc_poly_lj_module.f90 b/mc_poly_lj_module.f90 index eae764c..53217ae 100644 --- a/mc_poly_lj_module.f90 +++ b/mc_poly_lj_module.f90 @@ -38,22 +38,27 @@ MODULE mc_module INTEGER, PUBLIC :: n ! number of molecules REAL, DIMENSION(:,:), ALLOCATABLE, PUBLIC :: r ! centre of mass positions (3,n) REAL, DIMENSION(:,:), ALLOCATABLE, PUBLIC :: e ! quaternions (0:3,n) + REAL, DIMENSION(:,:,:), ALLOCATABLE, PUBLIC :: d ! bond vectors (0:3,na,n) - ! Private data - ! Bond vectors in body-fixed frame + ! Bond vectors in body-fixed frame (na and db are public) ! Isosceles triangle, 3 sites, with unit bond length and bond angle alpha, which we set to 75 degrees here - REAL, PARAMETER :: pi = 4.0*atan(1.0) - REAL, PARAMETER :: alpha = 75.0 * pi / 180.0, alpha2 = alpha / 2.0 - INTEGER, PARAMETER :: na = 3 - REAL, DIMENSION(3,na), PARAMETER :: db = RESHAPE ( [ & + REAL, PARAMETER :: pi = 4.0*ATAN(1.0) + REAL, PARAMETER :: alpha = 75.0 * pi / 180.0, alpha2 = alpha / 2.0 + INTEGER, PARAMETER, PUBLIC :: na = 3 + REAL, DIMENSION(3,na), PARAMETER, PUBLIC :: db = RESHAPE ( [ & & -SIN(alpha2), 0.0, -COS(alpha2)/3.0, & & 0.0, 0.0, 2.0*COS(alpha2)/3.0, & & SIN(alpha2), 0.0, -COS(alpha2)/3.0 ], [3,na] ) + ! Cutoff distance and force-shift parameters (all private) chosen as per the reference: + ! S Mossa, E La Nave, HE Stanley, C Donati, F Sciortino, P Tartaglia, Phys Rev E, 65, 041205 (2002) + REAL, PARAMETER :: r_cut = 2.612 ! in sigma=1 units, where r_cut = 1.2616 nm, sigma = 0.483 nm + REAL, PARAMETER :: sr_cut = 1.0/r_cut, sr_cut6 = sr_cut**6, sr_cut12 = sr_cut6**2 + REAL, PARAMETER :: lambda1 = 4.0*(7.0*sr_cut6-13.0*sr_cut12) + REAL, PARAMETER :: lambda2 = -24.0*(sr_cut6-2.0*sr_cut12)*sr_cut + INTEGER, PARAMETER :: lt = -1, gt = 1 ! Options for j-range - REAL, PARAMETER :: diameter = 2.0 * SQRT ( MAXVAL ( SUM(db**2,dim=1) ) ) ! Molecular diameter - ! Public derived type TYPE, PUBLIC :: potential_type ! A composite variable for interactions comprising REAL :: pot ! the potential energy cut and shifted to 0 at r_cut, and @@ -90,9 +95,10 @@ SUBROUTINE introduction IMPLICIT NONE INTEGER :: i + REAL :: diameter WRITE ( unit=output_unit, fmt='(a)' ) 'Lennard-Jones potential' - WRITE ( unit=output_unit, fmt='(a)' ) 'Cut-and-shifted' + WRITE ( unit=output_unit, fmt='(a)' ) 'Cut-and-force-shifted' WRITE ( unit=output_unit, fmt='(a)' ) 'Diameter, sigma = 1' WRITE ( unit=output_unit, fmt='(a)' ) 'Well depth, epsilon = 1' @@ -100,7 +106,12 @@ SUBROUTINE introduction DO i = 1, na ! Loop over atoms WRITE ( unit=output_unit, fmt='(a,i1,t40,3f15.6)' ) 'Body-fixed atom vector ', i, db(:,i) END DO ! End loop over atoms - WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Molecular diameter ', diameter + + diameter = 2.0 * SQRT ( MAXVAL ( SUM(db**2,dim=1) ) ) + WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Molecular diameter', diameter + WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'r_cut', r_cut + WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Force-shift lambda1', lambda1 + WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Force-shift lambda2', lambda2 END SUBROUTINE introduction @@ -111,15 +122,15 @@ SUBROUTINE conclusion END SUBROUTINE conclusion - SUBROUTINE allocate_arrays ( box, r_cut ) + SUBROUTINE allocate_arrays ( box ) IMPLICIT NONE REAL, INTENT(in) :: box ! simulation box length - REAL, INTENT(in) :: r_cut ! potential cutoff distance - REAL :: rm_cut_box + REAL :: rm_cut_box, diameter - ALLOCATE ( r(3,n), e(0:3,n) ) + ALLOCATE ( r(3,n), e(0:3,n), d(3,na,n) ) + diameter = 2.0 * SQRT ( MAXVAL ( SUM(db**2,dim=1) ) ) rm_cut_box = ( r_cut+diameter ) / box IF ( rm_cut_box > 0.5 ) THEN WRITE ( unit=error_unit, fmt='(a,f15.6)') 'rm_cut/box too large ', rm_cut_box @@ -131,15 +142,14 @@ END SUBROUTINE allocate_arrays SUBROUTINE deallocate_arrays IMPLICIT NONE - DEALLOCATE ( r, e ) + DEALLOCATE ( r, e, d ) END SUBROUTINE deallocate_arrays - FUNCTION potential ( box, r_cut ) RESULT ( total ) + FUNCTION potential ( box ) RESULT ( total ) IMPLICIT NONE TYPE(potential_type) :: total ! Returns a composite of pot, vir etc REAL, INTENT(in) :: box ! Simulation box length - REAL, INTENT(in) :: r_cut ! Potential cutoff distance ! total%pot is the nonbonded cut-and-shifted potential energy for whole system ! total%vir is the corresponding virial for whole system @@ -150,8 +160,13 @@ FUNCTION potential ( box, r_cut ) RESULT ( total ) TYPE(potential_type) :: partial ! Molecular contribution to total INTEGER :: i - IF ( SIZE(r,dim=2) /= n ) THEN ! should never happen - WRITE ( unit=error_unit, fmt='(a,2i15)' ) 'Array bounds error for r', n, SIZE(r,dim=2) + IF ( ANY ( SHAPE(r) /= [3,n] ) ) THEN ! should never happen + WRITE ( unit=error_unit, fmt='(a,3i15)' ) 'Array bounds error for r', n, SHAPE(r) + STOP 'Error in potential' + END IF + + IF ( ANY ( SHAPE(d) /= [3,na,n] ) ) THEN ! should never happen + WRITE ( unit=error_unit, fmt='(a,5i15)' ) 'Array bounds error for d', na, n, SHAPE(d) STOP 'Error in potential' END IF @@ -159,7 +174,7 @@ FUNCTION potential ( box, r_cut ) RESULT ( total ) DO i = 1, n - 1 - partial = potential_1 ( r(:,i), e(:,i), i, box, r_cut, gt ) + partial = potential_1 ( r(:,i), d(:,:,i), i, box, gt ) IF ( partial%ovr ) THEN total%ovr = .TRUE. ! Overlap detected @@ -174,16 +189,14 @@ FUNCTION potential ( box, r_cut ) RESULT ( total ) END FUNCTION potential - FUNCTION potential_1 ( ri, ei, i, box, r_cut, j_range ) RESULT ( partial ) - USE maths_module, ONLY : q_to_a + FUNCTION potential_1 ( ri, di, i, box, j_range ) RESULT ( partial ) IMPLICIT NONE - TYPE(potential_type) :: partial ! Returns a composite of pot, vir etc for given molecule - REAL, DIMENSION(3), INTENT(in) :: ri ! Coordinates of molecule of interest - REAL, DIMENSION(0:3), INTENT(in) :: ei ! Quaternion of molecule of interest - INTEGER, INTENT(in) :: i ! Index of molecule of interest - REAL, INTENT(in) :: box ! Simulation box length - REAL, INTENT(in) :: r_cut ! Potential cutoff distance - INTEGER, OPTIONAL, INTENT(in) :: j_range ! Optional partner index range + TYPE(potential_type) :: partial ! Returns a composite of pot, vir etc for given molecule + REAL, DIMENSION(3), INTENT(in) :: ri ! Coordinates of molecule of interest + REAL, DIMENSION(3,na), INTENT(in) :: di ! Bond vectors of molecule of interest + INTEGER, INTENT(in) :: i ! Index of molecule of interest + REAL, INTENT(in) :: box ! Simulation box length + INTEGER, OPTIONAL, INTENT(in) :: j_range ! Optional partner index range ! partial%pot is the nonbonded cut-and-shifted potential energy of molecule ri,ei with a set of other molecules ! partial%vir is the corresponding virial of molecule ri,ei @@ -194,26 +207,16 @@ FUNCTION potential_1 ( ri, ei, i, box, r_cut, j_range ) RESULT ( partial ) ! It is assumed that r has been divided by box ! Results are in LJ units where sigma = 1, epsilon = 1 - ! Note that this is the shifted LJ potential + ! Note that this is the force-shifted LJ potential with a linear smoothing term + ! S Mossa, E La Nave, HE Stanley, C Donati, F Sciortino, P Tartaglia, Phys Rev E, 65, 041205 (2002) INTEGER :: j, j1, j2, a, b - REAL :: rm_cut_box, rm_cut_box_sq, r_cut_sq - REAL :: sr2, sr6, sr12, rij_sq, rab_sq, virab, pot_cut + REAL :: diameter, rm_cut_box, rm_cut_box_sq, r_cut_sq + REAL :: sr2, sr6, sr12, rij_sq, rab_sq, virab, rmag REAL, DIMENSION(3) :: rij, rab, fab - REAL, DIMENSION(3,na) :: di, dj - REAL, DIMENSION(3,3) :: ai, aj REAL, PARAMETER :: sr2_ovr = 1.77 ! overlap threshold (pot > 100) TYPE(potential_type) :: pair - IF ( SIZE(r,dim=2) /= n ) THEN ! should never happen - WRITE ( unit=error_unit, fmt='(a,2i15)' ) 'Array bounds error for r', n, SIZE(r,dim=2) - STOP 'Error in potential_1' - END IF - IF ( SIZE(e,dim=2) /= n ) THEN ! should never happen - WRITE ( unit=error_unit, fmt='(a,2i15)' ) 'Array bounds error for e', n, SIZE(e,dim=2) - STOP 'Error in potential_1' - END IF - IF ( PRESENT ( j_range ) ) THEN SELECT CASE ( j_range ) CASE ( lt ) ! j < i @@ -231,24 +234,13 @@ FUNCTION potential_1 ( ri, ei, i, box, r_cut, j_range ) RESULT ( partial ) j2 = n END IF + diameter = 2.0 * SQRT ( MAXVAL ( SUM(db**2,dim=1) ) ) rm_cut_box = ( r_cut + diameter ) / box ! Molecular cutoff in box=1 units rm_cut_box_sq = rm_cut_box**2 ! squared r_cut_sq = r_cut**2 ! Potential cutoff squared in sigma=1 units - ! Calculate potential at cutoff - sr2 = 1.0 / r_cut_sq ! in sigma=1 units - sr6 = sr2**3 - sr12 = sr6**2 - pot_cut = sr12 - sr6 - partial = potential_type ( pot=0.0, vir=0.0, ovr=.FALSE. ) ! Initialize - ! Compute all space-fixed atom vectors for molecule i - ai = q_to_a ( ei ) ! Rotation matrix for i - DO a = 1, na ! Loop over all atoms - di(:,a) = MATMUL ( db(:,a), ai ) ! NB: equivalent to ai_T*db, ai_T=transpose of ai - END DO ! End loop over all atoms - DO j = j1, j2 ! Loop over selected range of partner molecules IF ( i == j ) CYCLE ! Skip self @@ -261,18 +253,12 @@ FUNCTION potential_1 ( ri, ei, i, box, r_cut, j_range ) RESULT ( partial ) rij = rij * box ! Now in sigma = 1 units - ! Compute all space-fixed atom vectors for molecule j - aj = q_to_a ( e(:,j) ) ! Rotation matrix for j - DO a = 1, na ! Loop over all atoms - dj(:,a) = MATMUL ( db(:,a), aj ) ! NB: equivalent to aj_T*db, aj_T=transpose of aj - END DO ! End loop over all atoms - ! Double loop over atoms on both molecules DO a = 1, na DO b = 1, na - rab = rij + di(:,a) - dj(:,b) ! Atom-atom vector, sigma=1 units - rab_sq = SUM ( rab**2 ) ! Squared atom-atom separation, sigma=1 units + rab = rij + di(:,a) - d(:,b,j) ! Atom-atom vector, sigma=1 units + rab_sq = SUM ( rab**2 ) ! Squared atom-atom separation, sigma=1 units IF ( rab_sq < r_cut_sq ) THEN ! Test within potential cutoff @@ -284,13 +270,13 @@ FUNCTION potential_1 ( ri, ei, i, box, r_cut, j_range ) RESULT ( partial ) RETURN ! Return immediately END IF + rmag = sqrt(rab_sq) sr6 = sr2**3 sr12 = sr6**2 - pair%pot = sr12 - sr6 ! LJ atom-atom pair potential (cut but not shifted) - virab = pair%pot + sr12 ! LJ atom-atom pair virial - pair%pot = pair%pot - pot_cut ! LJ atom-atom pair potential (cut-and-shifted) - fab = rab * virab * sr2 ! LJ atom-atom pair force - pair%vir = DOT_PRODUCT ( rij, fab ) ! Contribution to molecular virial + pair%pot = 4.0*(sr12-sr6) + lambda1 + lambda2*rmag ! LJ atom-atom pair potential (force-shifted) + virab = 24.0*(2.0*sr12-sr6 ) - lambda2*rmag ! LJ atom-atom pair virial + fab = rab * virab * sr2 ! LJ atom-atom pair force + pair%vir = DOT_PRODUCT ( rij, fab ) ! Contribution to molecular virial partial = partial + pair @@ -305,9 +291,8 @@ FUNCTION potential_1 ( ri, ei, i, box, r_cut, j_range ) RESULT ( partial ) END DO ! End loop over selected range of partner molecules ! Include numerical factors - partial%pot = partial%pot * 4.0 ! 4*epsilon - partial%vir = partial%vir * 24.0 / 3.0 ! 24*epsilon and divide virial by 3 - partial%ovr = .FALSE. ! No overlaps detected (redundant, but for clarity) + partial%vir = partial%vir / 3.0 ! Divide virial by 3 + partial%ovr = .FALSE. ! No overlaps detected (redundant, but for clarity) END FUNCTION potential_1