Merge branch 'master' of ssh://github.com/libAtoms/pymatnest

example_LJ_model.F90

@@ -55,18 +55,45 @@ subroutine ll_init_model(N_params, params)
    integer :: N_params
    double precision :: params(N_params)
 
-   epsilon(1,1) = 4.0
-   epsilon(2,2) = 4.0
-   epsilon(1,2) = 6.0
-   epsilon(2,1) = epsilon(1,2)
-
-   sigma(1,1) = 1.0
-   sigma(2,2) = 1.0
-   sigma(1,2) = (sigma(1,1)+sigma(2,2))/2.0
-   sigma(2,1) = sigma(1,2)
-
-   cutoff = 3.0*sigma
-   cutoff_sq = cutoff*cutoff
+   if (N_params==1) then
+
+       write(*,*) "No LJ parameters were found, the default will be used:"
+       epsilon(1,1) = 4.0
+       epsilon(2,2) = 4.0
+       epsilon(1,2) = 6.0
+       epsilon(2,1) = epsilon(1,2)
+
+       sigma(1,1) = 1.0
+       sigma(2,2) = 1.0
+       sigma(1,2) = (sigma(1,1)+sigma(2,2))/2.0
+       sigma(2,1) = sigma(1,2)
+
+       cutoff = 3.0*sigma
+       cutoff_sq = cutoff*cutoff
+
+   else
+
+       epsilon(1,1) = 4.0*params(1) ! The usual 4 factor in the LJ equation is
+                                    !  included here!
+       epsilon(2,2) = 4.0*params(2)
+       epsilon(1,2) = sqrt(epsilon(1,1)*epsilon(2,2))
+       epsilon(2,1) = epsilon(1,2)
+
+       sigma(1,1) = params(3)
+       sigma(2,2) = params(4)
+       sigma(1,2) = (sigma(1,1)+sigma(2,2))/2.0
+       sigma(2,1) = sigma(1,2)
+
+       cutoff = params(5)*sigma
+       cutoff_sq = cutoff*cutoff
+       write(*,*) "User defined LJ parameters:"
+       write(*,*) "(1)-(1)Epsilon=",params(1),"(2)-(2)Epsilon=",params(2)
+       write(*,*) "(1)-(1)Sigma=",  params(3),"(2)-(2)Sigma=",params(4)
+       write(*,*) "Mixed term LJ parameters are calculated using the", &
+                  &" Lorenz-Berthelot rule."
+       write(*,*) "LJ cutoff=",cutoff
+   endif
+
 
    E_offset = (sigma/cutoff)**12 - (sigma/cutoff)**6
 

fortranMCMDpy.py

@@ -95,6 +95,8 @@ def __init__(self, model_so):
            ctypes.c_void_p, # step_size_velo
            ctypes.c_void_p, # Emax
            ctypes.c_void_p, # nD
+           ctypes.c_void_p, # fixN
+           ctypes.c_void_p, # wall
            ctypes.c_void_p, # KEmax
            ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"), # final_E
            ndpointer(ctypes.c_int, flags="C_CONTIGUOUS"), # n_try
@@ -210,7 +212,7 @@ def MC_atom_walk_velo(self, at, n_steps, step_size, nD, KEmax):
         at.set_velocities(velo)
         return (n_try[0], n_accept[0], final_KE[0])
 
-    def MC_atom_walk(self, at, n_steps, step_size_pos, Emax, nD, KEmax=-1.0, step_size_velo=None):
+    def MC_atom_walk(self, at, n_steps, step_size_pos, Emax, nD, fixN=0, wall=0.0, KEmax=-1.0, step_size_velo=None):
         n = ctypes.c_int(len(at))
         n_steps = ctypes.c_int(n_steps)
         step_size_pos = ctypes.c_double(step_size_pos)
@@ -221,6 +223,8 @@ def MC_atom_walk(self, at, n_steps, step_size_pos, Emax, nD, KEmax=-1.0, step_si
         Emax = ctypes.c_double(Emax)
         KEmax = ctypes.c_double(KEmax)
         nD = ctypes.c_int(nD)
+        fixN = ctypes.c_int(fixN)
+        wall = ctypes.c_double(wall)
         pos = at.get_positions()
         if step_size_velo is None:
             velo = np.zeros( (1), dtype=np.float64 )
@@ -239,7 +243,7 @@ def MC_atom_walk(self, at, n_steps, step_size_pos, Emax, nD, KEmax=-1.0, step_si
         self.lib.fortran_mc_atom_(ctypes.byref(n), at.get_atomic_numbers().astype(np.int32), pos, velo, at.get_masses(), ctypes.byref(n_extra_data_c),
            extra_data, at.get_cell()[:,:],
            ctypes.byref(n_steps), ctypes.byref(step_size_pos), ctypes.byref(step_size_velo_c),
-           ctypes.byref(Emax), ctypes.byref(nD), ctypes.byref(KEmax), final_E, n_try, n_accept_pos, n_accept_velo)
+           ctypes.byref(Emax), ctypes.byref(nD), ctypes.byref(fixN), ctypes.byref(wall), ctypes.byref(KEmax), final_E, n_try, n_accept_pos, n_accept_velo)
         at.set_positions(pos)
         if n_extra_data_c.value > 0:
             at.arrays['ns_extra_data'][...] = extra_data

fortran_MC_MD.F90

@@ -64,20 +64,20 @@ subroutine fortran_MC_atom_velo(N, vel, mass, n_steps, step_size, nD, KEmax, fin
 end subroutine fortran_MC_atom_velo
 
 subroutine fortran_MC_atom(N, Z, pos, vel, mass, n_extra_data, extra_data, cell, n_steps, &
-                           step_size_pos, step_size_vel, Emax, nD, KEmax, final_E, n_try, n_accept_pos, n_accept_vel)
+                           step_size_pos, step_size_vel, Emax, nD, fixN, wall, KEmax, final_E, n_try, n_accept_pos, n_accept_vel)
    implicit none
    integer :: N
    integer :: Z(N)
    double precision :: pos(3,N), vel(3,N), mass(N), cell(3,3)
    integer :: n_extra_data
    double precision :: extra_data(n_extra_data,N)
    integer :: n_steps
-   double precision :: step_size_pos, step_size_vel, Emax, KEmax, final_E
-   integer :: n_try, n_accept_pos, n_accept_vel, nD
+   double precision :: step_size_pos, step_size_vel, Emax, KEmax, final_E, wall
+   integer :: n_try, n_accept_pos, n_accept_vel, nD, fixN
 
    logical :: do_vel
    integer :: d_i
-   double precision :: d_r, E, dE, KE, dKE, d_pos(3), d_vel(3)
+   double precision :: d_r, E, dE, KE, dKE, d_pos(3), d_vel(3), new_z
 
    double precision, external :: ll_eval_energy
    integer, external :: ll_move_atom_1
@@ -88,7 +88,7 @@ subroutine fortran_MC_atom(N, Z, pos, vel, mass, n_extra_data, extra_data, cell,
 
    do_vel = (step_size_vel /= 0.0)
 
-   n_try = N*n_steps
+   n_try = (N-fixN)*n_steps
    n_accept_pos = 0
    n_accept_vel = 0
    E = ll_eval_energy(N, Z, pos, n_extra_data, extra_data, cell)
@@ -102,7 +102,8 @@ subroutine fortran_MC_atom(N, Z, pos, vel, mass, n_extra_data, extra_data, cell,
       do i_at=1, N
          order(i_at) = i_at
       end do
-      do i_at=1, N-1
+      ! create random order in which single atom moves will be done
+      do i_at=1+fixN, N-1
          call random_number(d_r); d_i = floor(d_r*(N-i_at+1))+i_at
          if (d_i /= i_at) then
             t_i = order(i_at)
@@ -111,7 +112,8 @@ subroutine fortran_MC_atom(N, Z, pos, vel, mass, n_extra_data, extra_data, cell,
          endif
       end do
 
-      do i_at=1, N
+      ! go through list of moving atoms
+      do i_at=1+fixN, N
          d_i = order(i_at)
          if (do_vel) then
             call random_number(d_vel)
@@ -131,10 +133,27 @@ subroutine fortran_MC_atom(N, Z, pos, vel, mass, n_extra_data, extra_data, cell,
             endif
          endif
 
+         ! do single atom move of atom d_i
+         !    single atom move: only the change in energy contribution of the
+         !    single atom is calculated. 
          call random_number(d_pos)
          d_pos = 2.0*step_size_pos*(d_pos-0.5)
+
          if (nD==2) d_pos(3)=0.0
-         n_accept_pos = n_accept_pos + ll_move_atom_1(N, Z, pos, n_extra_data, extra_data, cell, d_i, d_pos, Emax-E, dE)
+
+         if (wall>0.0) then ! this is a surface simulation with a wall set
+            new_z=mod( (pos(3,d_i)+d_pos(3)), cell(3,3)) ! wrap Z coordinate
+            if (new_z < 0.0) new_z=new_z+cell(3,3) ! if coordinate was negative
+            if (new_z > (cell(3,3)-wall)) then
+                dE = 0.0 ! reject trial move, energy does not change
+            else ! proceed with the step
+                n_accept_pos = n_accept_pos + ll_move_atom_1(N, Z, pos, n_extra_data, extra_data, &
+                             & cell, d_i, d_pos, Emax-E, dE)
+            endif
+         else ! not a surface simulation, proceed as normal
+            n_accept_pos = n_accept_pos + ll_move_atom_1(N, Z, pos, n_extra_data, extra_data, cell, d_i, d_pos, Emax-E, dE)
+         endif
+
          E = E + dE
 
          if (do_vel .and.  vel_pos_rv >= 0.5) then
@@ -149,7 +168,9 @@ subroutine fortran_MC_atom(N, Z, pos, vel, mass, n_extra_data, extra_data, cell,
          endif
 
       end do
+
    end do
+   !write(*,*) "LIVIA", n_accept_pos, n_steps, step_size_pos, dE, E
 
    final_E = E
 

lammpslib.py

@@ -5,6 +5,7 @@
 import os
 import ctypes
 import operator
+import sys
 
 import numpy as np
 from numpy.linalg import norm
@@ -379,11 +380,11 @@ def set_lammps_pos(self, atoms):
 #        self.lmp.put_coosrds(lmp_c_positions)
         self.lmp.scatter_atoms('x', 1, 3, lmp_c_positions)
 
-
     def calculate(self, atoms, properties, system_changes):
         self.propagate(atoms, properties, system_changes, 0)
 
-    def propagate(self, atoms, properties, system_changes, n_steps, dt=None, dt_not_real_time=False, velocity_field=None):
+    def propagate(self, atoms, properties, system_changes, n_steps, dt=None,
+                  dt_not_real_time=False, velocity_field=None):
 
         """"atoms: Atoms object
             Contains positions, unit-cell, ...
@@ -404,15 +405,17 @@ def propagate(self, atoms, properties, system_changes, n_steps, dt=None, dt_not_
         if not self.started:
             self.start_lammps()
 
-        ####################################################################################################
-        #NB
+        ########################################################################
+        # NB
         if not self.initialized:
             self.initialise_lammps(atoms)
-        else: # still need to reset cell
-            # reset positions so that if they are cray from last propagation, change_box (in set_cell()) won't hang
-            # could do this only after testing for crazy positions?
-            # could also use scatter_atoms() to set values (requires MPI comm), or extra_atoms() to get pointers to local 
-            #     data structures to zero, but then will have to be careful with parallelism
+        else:  # Still need to reset cell
+            # Reset positions so that if they are crazy from last propagation,
+            # change_box (in set_cell()) won't hang.
+            # Could do this only after testing for crazy positions?
+            # Could also use scatter_atoms() to set values (requires MPI comm),
+            # or extra_atoms() to get pointers to local data structures to zero,
+            # but then will have to be careful with parallelism
             self.lmp.command("set atom * x 0.0 y 0.0 z 0.0")
             self.set_cell(atoms, change=True)
 
@@ -422,30 +425,32 @@ def propagate(self, atoms, properties, system_changes, n_steps, dt=None, dt_not_
         do_rebuild = False
         do_redo_atom_types = False
         try:
-            do_rebuild = ( len(atoms.numbers) != len(self.previous_atoms_numbers) ) or ( "numbers" in system_changes )
+            do_rebuild = (len(atoms.numbers) != len(self.previous_atoms_numbers)) or ("numbers" in system_changes)
             if not do_rebuild:
-                do_redo_atom_types = ( atoms.numbers != self.previous_atoms_numbers ).any()
+                do_redo_atom_types = (
+                        atoms.numbers != self.previous_atoms_numbers).any()
         except Exception:
            pass
 
-        self.lmp.command('echo none') # don't echo the atom positions
+        self.lmp.command('echo none')  # don't echo the atom positions
         if do_rebuild:
-           self.rebuild(atoms)
+            self.rebuild(atoms)
         elif do_redo_atom_types:
-           self.redo_atom_types(atoms)
-        self.lmp.command('echo log') # switch back log
+            self.redo_atom_types(atoms)
+        self.lmp.command('echo log')  # switch back log
 
         self.set_lammps_pos(atoms)
 
-        if n_steps > 0:
+        if n_steps > 0:  # TODO: here are velocities passed onto LAMMPS
             if velocity_field is None:
-                vel = atoms.get_velocities() / unit_convert("velocity", self.units)
+                vel = atoms.get_velocities() / unit_convert("velocity",
+                                                            self.units)
             else:
                 vel = atoms.arrays[velocity_field]
 
             # If necessary, transform the velocities to new coordinate system
             if self.coord_transform is not None:
-                vel = np.dot(self.coord_transform , np.matrix.transpose(vel) )
+                vel = np.dot(self.coord_transform, np.matrix.transpose(vel))
                 vel = np.matrix.transpose(vel)
 
             # Convert ase velocities matrix to lammps-style velocities array
@@ -454,9 +459,29 @@ def propagate(self, atoms, properties, system_changes, n_steps, dt=None, dt_not_
             # Convert that lammps-style array into a C object
             lmp_c_velocities =\
                 (ctypes.c_double * len(lmp_velocities))(*lmp_velocities)
-#            self.lmp.put_coosrds(lmp_c_velocities)
+            # self.lmp.put_coords(lmp_c_velocities)
             self.lmp.scatter_atoms('v', 1, 3, lmp_c_velocities)
 
+            # Keep atoms fixed
+            # # RY: use LAMMPS_init_cmds to set up NVE, 
+            # # e.g. group fixed id <= X; group mobile id > X; fix 1 mobile nve
+            # keep_atoms_fixed = int(sum([x == 0 for x in lmp_velocities]) / 3)
+            # if keep_atoms_fixed > 0:
+            #     self.lmp.command("group fixed id <= " + str(keep_atoms_fixed))
+            #     self.lmp.command("group mobile id > " + str(keep_atoms_fixed))
+                #self.lmp.command("fix freeze fixed setforce 0.0 0.0 0.0")
+                #if atoms.info["set_wall"]:
+                #    self.lmp.command("fix walls all wall/reflect zlo 0 zhi "
+                #                     + str(atoms.cell[2, 2]) + " units box")
+
+            # TODO: if we fix forces here, then it should be passed on, just
+            #  pass on keep_atoms_fixed
+            # TODO: if you have atoms with EXACTLY zero velocities, then freeze
+            #  them
+
+        # TODO: keep_atoms_fixed = 0 for potential energy calculations of the
+        #  initial configurations
+
         # Run for 0 time to calculate
         if dt is not None:
             if dt_not_real_time:
@@ -517,7 +542,7 @@ def propagate(self, atoms, properties, system_changes, n_steps, dt=None, dt_not_
         self.results['stress'] = stress * (-unit_convert("pressure", self.units))
 
 #        if 'forces' in properties:
-        f = np.zeros((len(atoms), 3))
+        f = np.zeros((len(atoms), 3))  # TODO: sets forces, doesn't update them
         f[:,:] = np.array([x for x in self.lmp.gather_atoms("f",1,3)]).reshape(-1,3)
         f *= unit_convert("force", self.units)
 
@@ -529,9 +554,11 @@ def propagate(self, atoms, properties, system_changes, n_steps, dt=None, dt_not_
         if not self.parameters.keep_alive:
             self.lmp.close()
 
-    def lammpsbc(self, pbc):
+    def lammpsbc(self, pbc, fix):
         if pbc:
             return 'p'
+        elif fix:
+            return 'f'
         else:
             return 's'
 
@@ -632,8 +659,13 @@ def initialise_lammps(self, atoms):
                 if 'boundary' in cmd:
                     break
             else:
-                self.lmp.command(
-                    'boundary ' + ' '.join([self.lammpsbc(bc) for bc in pbc]))
+                fix = False
+                # TODO: RBW – quick fix so that boundary parallel to surface
+                #  is not shrink wrapped
+                # if "set_wall" in atoms.info.keys():
+                #    fix = True
+                self.lmp.command('boundary ' + ' '.join([self.lammpsbc(bc, fix)
+                                                         for bc in pbc]))
 
         # Initialize cell
         self.set_cell(atoms, change=not self.parameters.create_box)
@@ -772,6 +804,7 @@ def initialise_lammps(self, atoms):
 
         self.initialized = True
 
+
 def write_lammps_data(filename, atoms, atom_types, comment=None, cutoff=None,
                       molecule_ids=None, charges=None, units='metal',
                       bond_types=None, angle_types=None, dihedral_types=None,

ns_analyse

@@ -92,6 +92,8 @@ def analyse(log_a, T, interval=1):
 
     Cvp = n_Extra_DOF*k_Boltzmann/2.0 + k_Boltzmann *beta*beta * (sum(Z_term * Es**2)/Z - U_pot**2)
 
+    # TODO: Modify script here to calculate ensemble-averaged order parameters
+
     if Vs is not None:
         V = sum(Z_term*Vs)/Z
         #thermal_exp = -1.0/V * (sum(Z_term*Vs*Vs)*(-beta)*Z - sum(Z_term*Vs)*sum(Z_term*Vs)*(-beta)) / Z**2