!!!!! This program calculate heat capacity for water at a specific temperature !!!!! !!!!! User have to modify values in MODIFY SPACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Program heatcapacity double precision,dimension(3000000)::t,temp,vol,press,ke,pe,e double precision:: em, ems, var, cv, kb, na, TK, mol, temper, volume, eps, sig, cal, pressure, omass, hmass, wmass, density integer::n,no !!! constants: kb=1.3806503*(10.**(-23.)) na=6.02214199*(10.**(23.)) cal=4.18400 !!!!! MODIFY SPACE !!!!! n=10000 ! total time steps no=2001 ! time steps where the system starts equilibrate omass=15.9994 ! mass of an oxygen atom hmass=1.008 ! mass of a hydorgen atom mol=1000 ! no. of water molecules open(1,file='tip4p.nvt4.out') !!!!!!!!!!!!!!!!!!!!!!!! do i=1, n read(1,*) t(i),temp(i),vol(i),press(i),ke(i),pe(i),e(i) end do close(1) !!! calculate heatcapacity em=0. ems=0. temper=0. volume=0. pressure=0. wmass=((1.*omass+2.*hmass)/na)/1000. do i=no, n !convert to SI unit e(i)=e(i)*cal*1000./na vol(i)=vol(i)*((1.0*(10.**(-10.)))**(3.)) pressure=pressure+press(i) temper=temper+temp(i) volume=volume+vol(i) em=em+e(i) ems=ems+(e(i)**2.) end do pressure=pressure/((n-no+1)*1.) temper=temper/((n-no+1)*1.) volume=volume/((n-no+1)*1.) var= ems/((n-no+1)*1.) - (em/((n-no+1)*1.))**2. cv=var/(kb*(temper**2.)) density=mol*wmass/volume !convert to KJ K^-1 kg^-1 cv=(cv/(mol*wmass))/1000.0 write(*,*) 'cv = ', cv write(*,*) 'temperature = ', temper write(*,*) 'density = ', density write(*,*) 'pressure = ', pressure stop end program heatcapacity