Differences

This shows you the differences between two versions of the page.

--- exercises:2018_uzh_acpc2:mol_sol [2018/05/03 13:00] – jglan
+++ exercises:2018_uzh_acpc2:mol_sol [2020/08/21 10:15] (current) – external edit 127.0.0.1
@@ Line 4: / Line 4: @@
 [[http://www1.lsbu.ac.uk/water/water_models.html|Water molecular models]] are computational techniques that have been developed in order to help discover the structure of water. In this section, you will be asked to calculate some physical properties based on classical molecular dynamics simulation. The TIP3/Fw model will be usded in the simulations.
-We have prepared a CP2K input file ''water.inp'' for running a MD simulation of liquid water using the force field from the first exercise (parametrized by [[https://aip.scitation.org/doi/pdf/10.1063/1.1884609|Praprotnik et al.]]). Download {{ water.zip |water.zip}}
+We have prepared a CP2K input file ''water.inp'' for running a MD simulation of liquid water using the force field from the first exercise (parametrized by [[https://aip.scitation.org/doi/pdf/10.1063/1.1884609|Praprotnik et al.]]).
+<note important>Download {{ water.zip |water.zip}} and extract it.</note>
@@ Line 32: / Line 33: @@
   * Plot RMSD for the water at 300K and calculate corresponding diffusion coefficient, are they expected?
 </note>
+<note important>The diffusion coefficient is calculated using MSD but NOT RMSD.</note>
 We will compute the vibrational spectrum, and dielectric constant of water based on molecular dynamics. The spectra for water are available in this paper [[https://aip.scitation.org/doi/pdf/10.1063/1.1884609|https://doi.org/10.1063/1.1884609]]. The provided program computes the correlation function of the (derivative of) the dipole moment and performs the Fourier transform.
@@ Line 53: / Line 55: @@
 \operatorname{Var}(M) = (\langle M \cdot M\rangle - \langle M\rangle\langle M\rangle ) \ .
 \end{equation}
+Compile the FORTRAN code, and execute the program
+<code>
+gfortran cpt_ir_diele.f90 -o cpt_ir_diele.o
+./cpt_ir_diele.o < dipole.in
+</code>
 <note>**TASK 4**
@@ Line 60: / Line 68: @@
 </note>
-<code>
-!
-!
-IMPLICIT NONE
-INTEGER, PARAMETER :: dp=KIND(0.0D0)
-REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: correlation
-REAL(KIND=dp) :: integral,omega,Pi,timestep
-REAL(KIND=dp) :: M,M_aver,VAR,epsilon_0,e,kb,angstrom2meter,s,epsilon,v,t
-REAL(KIND=dp), DIMENSION(3,1000000) :: dipder,dip
-REAL(KIND=dp), DIMENSION(3) :: m_vec
-INTEGER :: N,I,J,Nmax,Fmax
-CHARACTER(LEN=100) :: line,filename
-Pi=4.0D0*ATAN(1.0D0)
-READ(5,*) filename
-READ(5,*) timestep
-OPEN(10,FILE=filename)
-N=0
-DO
- READ(10,'(A100)',END=999) line
- IF (INDEX(line,' MM DIPOLE [BERRY PHASE](Debye)|').NE.0) THEN
-   ! N=N+1
-    READ(line(45:),*) dip(:,N)
- ENDIF
- IF (INDEX(line,' MM DIPOLE [BERRY PHASE] DERIVATIVE(A.U.)|').NE.0) THEN
-    N=N+1
-    READ(line(45:),*) dipder(:,N)
- ENDIF
-ENDDO
-CONTINUE
-CLOSE(10)
-! use only the first 10% of the data for the correlation function, as the rest is not statistically meaningful
-! and our quadratic algorithm becomes too slow
-Nmax=N/10
-print *, Nmax
-ALLOCATE(correlation(0:Nmax))
-correlation=0.0_dp
-DO I=1,N-Nmax
- DO J=I,I+Nmax
-    correlation(J-I)=correlation(J-I)+DOT_PRODUCT(dipder(:,I),dipder(:,J))
- ENDDO
-ENDDO
-DO I=0,Nmax
-   correlation(I)=correlation(I)/(REAL(N-I,kind=dp)*REAL(N,kind=dp))
-ENDDO
-OPEN(UNIT=10,FILE="dip_dip_correlation.time")
-write(10,*) "# correlation in the time domain, first column time in fs"
-DO I=-Nmax,Nmax
-   write(10,*) I*timestep,correlation(ABS(I))/correlation(0)
-ENDDO
-CLOSE(10)
-OPEN(UNIT=10,FILE="dip_dip_correlation.freq")
-write(10,*) "# correlation in the frequency domain [cm^-1]"
-!Fmax up to 4000 cm^-1
-Fmax=N*((4000D0/(2*Pi))*(2.0D0*Pi*timestep*1.0D-15*29979245800.0))
-DO I=0,Fmax
-   omega=(2.0D0*Pi*I)/N
-   integral=0.0_dp
-   DO J=0,Nmax
-      integral=integral+cos(omega*J)*correlation(J)
-   ENDDO
-   ! frequency in cm^-1
-   write(10,*) omega/(2.0D0*Pi*timestep*1.0D-15*29979245800.0),integral
-ENDDO
-CLOSE(10)
-M=0
-M_aver=0
-dip = dip * 3.33564 * 1e-30  ! convert Deybe to C*m [SI]
-DO I=1,N
-  M = M + DOT_PRODUCT(dip(:,I),dip(:,I))
-  M_vec(:) = M_vec(:) + dip(:,I)
-ENDDO
-  M_aver = sqrt(DOT_PRODUCT(M_vec(:),M_vec(:)))
-  M_aver = M_aver*M_aver / N
-  VAR = ( M - M_aver ) / (N-1)
-!write(*,*) var
-T=300
-epsilon_0 = 8.8541878176e-12 ![F/m] vacuum permittivity
-e = 1.602176565e-19 ! [C] elementary charge
-kb = 1.3806488e-23 ![J/K] Boltzmann constant
-angstrom2meter = 1e-30 ! volume
-V=9.865*9.865*9.865
-s = (4*pi)/(3*V*kb*T*angstrom2meter*epsilon_0)
-epsilon = 1 + s*var
-!write(*,*) "M,M_aver",M, M_aver,N
-write(*,*) "Dielectric Constant",epsilon !,s,var
-END
-</code>
 ===== Ramachandran plot =====
@@ Line 184: / Line 85: @@
 Visualize the structure ''glyala.pdb'' with VMD and determine the atomic indices of the atoms defining the dihedral angles.
 </note>
+<note important>//Note:// While VMD starts counting atoms from 0, CP2K starts counting from 1, i.e. the VMD indices need to be increased by 1.</note>
 With this knowledge at hand,
@@ Line 190: / Line 94: @@
 <note>**TASK 2**
-  - The atomic indices defining the dihedral indices in the input file ''geo.in'' are missing. Replace ''I1'' to ''I4'' by the atomic indices determined previously. //Note:// While VMD starts counting atoms from 0, CP2K starts counting from 1, i.e. the VMD indices need to be increased by 1.
+  - The atomic indices defining the dihedral indices in the input file ''geo.in'' are missing. Replace ''I1'' to ''I4'' by the atomic indices determined previously.
   - Use ''perform-gopt.sh'' to perform the grid of geometry optimizations.
   - Use gnuplot to plot the potential energy surface (we have provided a script ''epot.gp''). Which are the two most favoured conformations? <code> $ gnuplot</code><code> gnuplot > load "epot.gp"</code>
 </note>
 ===== Glyala in water =====
 Now, we will move to a more realistic system - Glyala in water. We will preformed a MD of glyala in water and save the trajectory.
@@ Line 203: / Line 106: @@
 <note>**TASK 6**
-   - Perform the molecular dynamics simulation using NVT ensemble at 300K.
+   - Perform the molecular dynamics simulation using NVT ensemble at 300K. Change TIMECON (i.e.500, 2000 fs) in the &THERMOSTAT section.
-   - Re-run the calculation using NVT ensemble with different TIMECON (500, 2000 fs) in the &THERMOSTAT section, and plot the total energy, temperature against time. Explain what you observe.
    - Determine from which step the system is equilibrated, plot the calculated properties and explain why.
    - Compute the O-O radial distribution function for water with acceptable statistics using 20 ps (after equilibration) of simulated time.
+   - Determine the solvation shell by calculating RDF of g$_{CO}$ (carbon atoms from glyala and oxygen atoms from water)
 </note>