exercises:2018_uzh_acpc2:mol_sol
Differences
This shows you the differences between two versions of the page.
Both sides previous revisionPrevious revisionNext revision | Previous revision | ||
exercises:2018_uzh_acpc2:mol_sol [2018/05/03 12:53] – jglan | exercises:2018_uzh_acpc2:mol_sol [2020/08/21 10:15] (current) – external edit 127.0.0.1 | ||
---|---|---|---|
Line 4: | Line 4: | ||
[[http:// | [[http:// | ||
- | We have prepared a CP2K input file '' | + | We have prepared a CP2K input file '' |
+ | <note important> | ||
Line 32: | Line 33: | ||
* Plot RMSD for the water at 300K and calculate corresponding diffusion coefficient, | * Plot RMSD for the water at 300K and calculate corresponding diffusion coefficient, | ||
</ | </ | ||
+ | <note important> | ||
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:// | 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:// | ||
Line 53: | Line 55: | ||
\operatorname{Var}(M) = (\langle M \cdot M\rangle - \langle M\rangle\langle M\rangle ) \ . | \operatorname{Var}(M) = (\langle M \cdot M\rangle - \langle M\rangle\langle M\rangle ) \ . | ||
\end{equation} | \end{equation} | ||
+ | |||
+ | Compile the FORTRAN code, and execute the program | ||
+ | < | ||
+ | gfortran cpt_ir_diele.f90 -o cpt_ir_diele.o | ||
+ | ./ | ||
+ | </ | ||
< | < | ||
Line 60: | Line 68: | ||
</ | </ | ||
- | < | ||
- | ! | ||
- | ! | ||
- | IMPLICIT NONE | ||
- | INTEGER, PARAMETER :: dp=KIND(0.0D0) | ||
- | REAL(KIND=dp), | ||
- | REAL(KIND=dp) :: integral, | ||
- | REAL(KIND=dp) :: M, | ||
- | REAL(KIND=dp), | ||
- | REAL(KIND=dp), | ||
- | INTEGER :: N, | ||
- | CHARACTER(LEN=100) :: line, | ||
- | |||
- | Pi=4.0D0*ATAN(1.0D0) | ||
- | |||
- | READ(5,*) filename | ||
- | READ(5,*) timestep | ||
- | |||
- | OPEN(10, | ||
- | N=0 | ||
- | DO | ||
- | | ||
- | IF (INDEX(line,' | ||
- | ! N=N+1 | ||
- | READ(line(45: | ||
- | ENDIF | ||
- | IF (INDEX(line,' | ||
- | N=N+1 | ||
- | READ(line(45: | ||
- | ENDIF | ||
- | ENDDO | ||
- | 999 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: | ||
- | correlation=0.0_dp | ||
- | |||
- | DO I=1,N-Nmax | ||
- | DO J=I,I+Nmax | ||
- | correlation(J-I)=correlation(J-I)+DOT_PRODUCT(dipder(:, | ||
- | ENDDO | ||
- | ENDDO | ||
- | |||
- | DO I=0,Nmax | ||
- | | ||
- | ENDDO | ||
- | |||
- | OPEN(UNIT=10, | ||
- | write(10,*) "# correlation in the time domain, first column time in fs" | ||
- | DO I=-Nmax, | ||
- | | ||
- | ENDDO | ||
- | CLOSE(10) | ||
- | |||
- | OPEN(UNIT=10, | ||
- | write(10,*) "# correlation in the frequency domain [cm^-1]" | ||
- | !Fmax up to 4000 cm^-1 | ||
- | Fmax=N*((4000D0/ | ||
- | DO I=0,Fmax | ||
- | | ||
- | | ||
- | DO J=0,Nmax | ||
- | integral=integral+cos(omega*J)*correlation(J) | ||
- | ENDDO | ||
- | ! frequency in cm^-1 | ||
- | | ||
- | 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(:, | ||
- | M_vec(:) = M_vec(:) + dip(:,I) | ||
- | ENDDO | ||
- | M_aver = sqrt(DOT_PRODUCT(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)/ | ||
- | epsilon = 1 + s*var | ||
- | |||
- | !write(*,*) " | ||
- | write(*,*) " | ||
- | |||
- | |||
- | END | ||
- | </ | ||
===== Ramachandran plot ===== | ===== Ramachandran plot ===== | ||
Line 184: | Line 85: | ||
Visualize the structure '' | Visualize the structure '' | ||
</ | </ | ||
+ | |||
+ | <note important>// | ||
+ | |||
With this knowledge at hand, | With this knowledge at hand, | ||
Line 190: | Line 94: | ||
< | < | ||
- | - The atomic indices defining the dihedral indices in the input file '' | + | - The atomic indices defining the dihedral indices in the input file '' |
- Use '' | - Use '' | ||
- Use gnuplot to plot the potential energy surface (we have provided a script '' | - Use gnuplot to plot the potential energy surface (we have provided a script '' | ||
</ | </ | ||
- | |||
===== Glyala in water ===== | ===== 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. | 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: | ||
< | < | ||
- | - Perform the molecular dynamics simulation using NVT ensemble at 300K. | + | - Perform the molecular dynamics simulation using NVT ensemble at 300K. Change |
- | - Re-run the calculation using NVT ensemble with different | + | |
- Determine from which step the system is equilibrated, | - Determine from which step the system is equilibrated, | ||
- Compute the O-O radial distribution function for water with acceptable statistics using 20 ps (after equilibration) of simulated time. | - 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) | ||
</ | </ | ||
exercises/2018_uzh_acpc2/mol_sol.1525351980.txt.gz · Last modified: 2020/08/21 10:15 (external edit)