xiaoquan

MD.log

@@ -0,0 +1,14 @@
+
+                            < M A T L A B (R) >
+                  Copyright 1984-2015 The MathWorks, Inc.
+                   R2015a (8.5.0.197613) 64-bit (glnxa64)
+                             February 12, 2015
+
+
+For online documentation, see http://www.mathworks.com/support
+For product information, visit www.mathworks.com.
+
+
+	Academic License
+
+>> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> 

MD.m

@@ -0,0 +1,148 @@
+clear all;
+
+%Build output files.
+%trajectory.xyz is the trajectory file without PBC(Periodic Boundary Condtion).
+%trajectory_imaged.xyz is the trajectory file with PBC.
+%md.out is the output file contains information of Step, Temperature, Energies(Kinetic, Potential, Total).
+%msd_r.out is MSD(Mean Square Displacement) for diffusion coefficient calculation.
+s=warning;
+warning off all;
+delete('trajectory.xyz','trajectory_imaged.xyz','md.out','msd_r.out');
+warning(s);
+file_trajectory=fopen('trajectory.xyz','a');
+file_trajectory_imaged=fopen('trajectory_imaged.xyz','a');
+file_energy=fopen('md.out','a');
+file_msd=fopen('msd_r.out','a');
+
+%Parameters.
+kb=1;%Boltzmann constant.
+natom=100;%Number of atoms.
+boxsize=14;%The side-length of the cubic box, unit angstrom.
+mass=1;%Mass of one atom.
+dt=1e-3;%Time step. Unit second.
+temp=300;%Temperature, unit K.
+step=100000;%Total running step.
+outstep=50;%Output frequency.
+epsilon=1;%Parameter for LJ potential, the depth of the potential well.
+sigma=2.5;%Parameter for LJ potential, the distance at which the inter-particle potential is 0.
+rc=2.5*sigma;%Cutoff distance, the LJ potential is truncated at the cutoff distance.
+ecut=4*epsilon*((sigma/rc)^12 - (sigma/rc)^6);%LJ potential at the cutoff distance.
+
+tic; %Timer
+%Initilization position.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+n=floor(power(natom,1/3))+1;
+for i=1:power(n,3)
+    p(1,i)=mod(i,n);
+    if (mod(i,n)==0)
+        p(1,i)=n;
+    end
+end
+    j=1;
+for i=1:n:power(n,3)
+    j=mod(j,n);
+    if (mod(j,n)==0)
+        j=n;
+    end
+    p(2,i:i+(n-1))=j;
+    j=j+1;
+end
+    j=1;
+for i=1:power(n,2):power(n,3)
+    p(3,i:i+(power(n,2)-1))=j;
+    j=j+1;
+end
+p=p*boxsize/n;%Adjust the lattice to the box.
+p=p(:,randperm(power(n,3),natom));%p is the position for all atoms.
+p_ref=p; %Choose initial position as reference for diffusion coefficient calculation.
+
+%Calculate potential and force for the initial condition.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+pen=0;%pen is the potential energy.
+f=zeros(3,natom);%f is the force.
+for i=1:(natom-1)
+    for j=(i+1):natom
+        r1=p(:,i)-p(:,j);
+        r1=r1-boxsize*round(r1./boxsize); %Consider PBC to get the correct distance between two atoms.
+        r=sqrt(sum(r1.^2));
+        if (r <= rc)
+            ff=24*epsilon*sigma^6*(2*sigma^6-r^6)/(r^14);
+            f(:,i)=f(:,i)+ff*r1;
+            f(:,j)=f(:,j)-ff*r1;
+            pen=pen+4*epsilon*((sigma/r)^12-(sigma/r)^6)-ecut;
+        end
+    end
+end
+acce=f/mass;%Calculate acceleration on each atom.
+
+%Initilization velocity.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+v=rand(3,natom); %Generate random numbers as initial velocities.
+ave_v=sum(v,2)/natom; %Calculate the average value of velocities.
+v=(v-repmat(ave_v,[1,natom])); %Scale velocities to make sure the total momentum is zero.
+ave_v=sum(v,2)/natom;
+ave_v2=sum(sum(v.^2,2)/natom);
+fa=sqrt(3*kb*temp/(mass*ave_v2)); %Velocity scale factor to specific temperature.
+v=(v-repmat(ave_v,[1,natom]))*fa; %Scale velocites to make sure the total momentum is zero and fit the temperature.
+ken=sum(sum(0.5*mass*v.^2,2)); %ken is the kinetic energy.
+
+%Integration.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+for m=1:step
+    p=p+v*dt+0.5*acce*dt^2; %Update position.
+
+    %Calculate force and potential according to the new position.
+    pen=0;
+    f=zeros(3,natom);
+    for i=1:(natom-1)
+        for j=(i+1):natom
+            r1=p(:,i)-p(:,j);
+            r1=r1-boxsize*round(r1./boxsize);
+            r=sqrt(sum(r1.^2));
+            if (r <= rc)
+                ff=24*epsilon*sigma^6*(2*sigma^6-r^6)/(r^14);
+                f(:,i)=f(:,i)+ff*r1;
+                f(:,j)=f(:,j)-ff*r1;
+                pen=pen+4*epsilon*((sigma./r).^12-(sigma./r).^6)-ecut;
+            end
+        end
+    end
+
+    v=v+0.5*acce*dt+0.5*f/mass*dt^2;%Update velocity.
+    ken=sum(sum(0.5*mass*v.^2,2));%Calculate kinetic energy.
+    realt=mass*sum(sum(v.^2,2)/natom)/(3*kb);%Calculate real temperature.
+
+    %Output according to the frequency setting.
+    if (mod(m,outstep)==0) %When the running step fit the output frequency, output the following information.
+        %Output MSD result.
+        msd_r=sum(sum((p-p_ref).^2,2))/natom; %calculte MSD
+        fprintf(file_msd,[num2str(m),'\t',num2str(msd_r,'%5.3e'),'\n']);
+
+        %Output trajectory file without PBC.
+        fprintf(file_trajectory,[num2str(natom),'\n']);
+        fprintf(file_trajectory,[num2str(m),' step \n']);
+        for i = 1:natom
+            fprintf(file_trajectory,['C','\t']);
+            for j =1:3
+                fprintf(file_trajectory,[num2str(p(j,i)'),'\t']);
+            end
+            fprintf(file_trajectory,'\n');
+        end
+
+        %Output trajectory file with PBC.
+        pi=p-boxsize*round(p./boxsize); %pi is the postion with PBC.
+        fprintf(file_trajectory_imaged,[num2str(natom),'\n']);
+        fprintf(file_trajectory_imaged,[num2str(m),' step \n']);
+        for i = 1:natom
+            fprintf(file_trajectory_imaged,['C','\t']);
+            for j =1:3
+                fprintf(file_trajectory_imaged,[num2str(pi(j,i)'),'\t']);
+            end
+            fprintf(file_trajectory_imaged,'\n');
+        end
+
+        %Output step and energies.
+        fprintf(file_energy,[num2str(m),'\t','T= ',num2str(realt,'%5.2f'),'\t']);
+        fprintf(file_energy,['Kinetic= ',num2str(ken,'%5.3e'),'\t']);
+        fprintf(file_energy,['Potential= ',num2str(pen,'%5.3e'),'\t']);
+        fprintf(file_energy,['Total= ',num2str((ken+pen),'%5.3e'),'\n']);
+    end
+end
+fprintf(file_energy,['#', num2str(step),' steps MD simulations cost ',num2str(toc),' seconds. \n']); %Output the total running time
+fprintf(file_energy,['#', 'The speed is about ',num2str(step/toc,'%8.3f'),' step/second. \n']); %Output the calculation speed.