MATLAB Examples

# solvate_atom.m

## Function arguments

• [limits] is a 1x3 or 1x6 volume variable
• density (number) is the density
• r is a number and the minimum distance between solvent-solute particles.
• maxsol (number) maxsol is the max number of solvent molecules
• solute_atom is the existing atom struct which needs solvation
• Optional string can be the desired water model, like 'SPC' or 'TIP4P', or with 'custom' be used to solvate with a custom structure, like a ethanol slab

## Examples

• SOL = solvate_atom(limits,density,r,maxsol,solute_atom)
• SOL = solvate_atom(limits,density,r,maxsol,solute_atom,'tip4p')
• SOL = solvate_atom(limits,density,r,maxsol,solute_atom,'spc_ice')

## Dependencies

• import_atom_gro
• scale_atom
• replicate_atom
• translate_atom
• merge_atom
• slice_atom
• update_atom
```function SOL = solvate_atom(limits,density,r,maxsol,solute_atom,varargin) if numel(limits)==1 limits(4:6)=limits(1); limits(1:3)=0; Lx=limits(4); Ly=limits(5); Lz=limits(6); elseif numel(limits)==3 Lx=limits(1); Ly=limits(2); Lz=limits(3); limits(4)=limits(1); limits(5)=limits(2); limits(6)=limits(3); limits(1:3)=0; elseif numel(limits)==6 Lx=limits(4)-limits(1); Ly=limits(5)-limits(2); Lz=limits(6)-limits(3); end % Old way of doing it... % SOL=import_atom_gro(strcat('spc864_',num2str(density,'%3.2f'),'.gro')); % I was fastest... [2 2 1]*216.gro if nargin>5 watermodel=varargin(1); if iscell(watermodel) assignin('caller','watermodel',watermodel) watermodel=char(watermodel{1}); end if strncmpi(watermodel,'spc_ice',7) SOL=import_atom_gro('96spc_hex_ice_h.gro'); pause(1) disp('Adding hexagonal spc ice!!!') elseif strncmpi(watermodel,'tip4p_ice',9) SOL=import_atom_gro('96tip4p_hex_ice_h.gro'); pause(1) disp('Adding hexagonal tip4p ice!!!') elseif strncmpi(watermodel,'spce',4) SOL=import_atom_gro('864_spce.gro'); elseif strncmpi(watermodel,'spc',3) SOL=import_atom_gro('864_spc.gro'); elseif strncmpi(watermodel,'tip3p',5) SOL=import_atom_gro('864_tip3p.gro'); elseif strncmpi(watermodel,'tip4p',5) SOL=import_atom_gro('864_tip4p.gro'); elseif strncmpi(watermodel,'tip5p',5) SOL=import_atom_gro('864_tip5p.gro'); elseif strncmpi(watermodel,'custom',5) SOL=varargin{2}; if nargin>7 Box_dim=varargin{3}; end end else SOL=import_atom_gro('864_spc.gro'); end atomsperSOL=sum([SOL.molid]==1); SOL=scale_atom(SOL,Box_dim,[1 1 1]./density,'all') nx=ceil(Lx/Box_dim(1)) ny=ceil(Ly/Box_dim(2)) nz=ceil(Lz/Box_dim(3)) SOL=replicate_atom(SOL,Box_dim,[nx ny nz],'xzy'); disp('Replicated n times'); nx*ny*nz; if (limits(1)+limits(2)+limits(3)) ~= 0 disp('Translating the water box'); SOL=translate_atom(SOL,[limits(1) limits(2) limits(3)],'all'); end % assignin('caller','SOL',SOL); % pause disp('nSOL before merge'); size(SOL,2)/atomsperSOL; if size(solute_atom,2) > 0 if size(SOL,2) > 10000 || size(solute_atom,2) > 10000 nSOL_block=size(SOL,2)/(nx*ny*nz); SOL_count=1;SOL_merged=[];count=1; while SOL_count< size(SOL,2) SOL_block= SOL(SOL_count:SOL_count+nSOL_block-1); SOL_block = merge_atom(solute_atom,limits(4:6)-.4,SOL_block,'molid','Hw',[r-.4 r]); SOL_merged = [SOL_merged SOL_block]; SOL_count=SOL_count+nSOL_block; disp('Number of water molecules...'); count=count+1; size(SOL_merged,2)/atomsperSOL end SOL=SOL_merged; else SOL = merge_atom(solute_atom,limits(4:6)-.4,SOL,'molid','Hw',[r-.4 r]); end else SOL = slice_atom(SOL,[limits(1) limits(2) limits(3) limits(4:6)-.4],0); end SOL=update_atom(SOL); % Randomize the order of the SOL molecules nSOL=size(SOL,2); if iscellstr({maxsol}) == 1 maxsol=nSOL/atomsperSOL; end rand_molid=randperm(nSOL/atomsperSOL); rand_molid=rand_molid(1:maxsol); rand_index=ismember([SOL.molid],rand_molid); SOL=SOL(rand_index); % Delete water molecules if not using the <maxsol> option if iscellstr({maxsol}) == 0 if atomsperSOL*maxsol > size(SOL,2) disp('Ooops, you asked for too much solvent...') disp('Maximum number of solvent molecules allowed without changing the density or rmin is:') size(SOL,2)/atomsperSOL SOL=SOL(1:atomsperSOL*maxsol); else % Use randperm instead... ? SOL=SOL(1:atomsperSOL*maxsol); end end SOL=update_atom(SOL); disp('nSOL after merge') size(SOL,2)/atomsperSOL assignin('caller','SOL',SOL); assignin('caller','limits',limits); ```