# 1 "f.c" # 1 "/usr/include/stdio.h" 1 3 # 1 "/usr/include/libio.h" 1 3 # 1 "/usr/include/features.h" 1 3 # 117 "/usr/include/features.h" 3 # 1 "/usr/include/sys/cdefs.h" 1 3 # 1 "/usr/include/features.h" 1 3 # 222 "/usr/include/features.h" 3 # 22 "/usr/include/sys/cdefs.h" 2 3 # 52 "/usr/include/sys/cdefs.h" 3 # 86 "/usr/include/sys/cdefs.h" 3 # 205 "/usr/include/features.h" 2 3 # 1 "/usr/include/gnu/stubs.h" 1 3 # 219 "/usr/include/features.h" 2 3 # 29 "/usr/include/libio.h" 2 3 # 1 "/usr/include/_G_config.h" 1 3 # 1 "/usr/include/gnu/types.h" 1 3 typedef unsigned char __u_char; typedef unsigned short __u_short; typedef unsigned int __u_int; typedef unsigned long __u_long; typedef struct { long int __val[2]; } __quad_t; typedef struct { __u_long __val[2]; } __u_quad_t; typedef __quad_t *__qaddr_t; typedef __u_quad_t __dev_t; typedef __u_int __uid_t; typedef __u_int __gid_t; typedef __u_long __ino_t; typedef __u_int __mode_t; typedef __u_int __nlink_t; typedef long int __off_t; typedef __quad_t __loff_t; typedef int __pid_t; typedef int __ssize_t; typedef struct { int __val[2]; } __fsid_t; typedef int __daddr_t; typedef char *__caddr_t; typedef long int __time_t; typedef long int __swblk_t; typedef long int __clock_t; typedef unsigned long int __fd_mask; typedef struct { __fd_mask fds_bits[1024 / (8 * sizeof (__fd_mask)) ]; } __fd_set; typedef int __key_t; typedef unsigned short int __ipc_pid_t; # 9 "/usr/include/_G_config.h" 2 3 # 1 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 1 3 # 19 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 61 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 131 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 typedef unsigned int size_t; # 258 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 typedef unsigned int wint_t; # 302 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 13 "/usr/include/_G_config.h" 2 3 typedef int _G_int16_t ; typedef int _G_int32_t ; typedef unsigned int _G_uint16_t ; typedef unsigned int _G_uint32_t ; # 31 "/usr/include/libio.h" 2 3 # 1 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stdarg.h" 1 3 typedef void *__gnuc_va_list; # 116 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stdarg.h" 3 # 202 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stdarg.h" 3 # 47 "/usr/include/libio.h" 2 3 # 64 "/usr/include/libio.h" 3 # 98 "/usr/include/libio.h" 3 struct _IO_jump_t; struct _IO_FILE; typedef void _IO_lock_t; struct _IO_marker { struct _IO_marker *_next; struct _IO_FILE *_sbuf; int _pos; # 182 "/usr/include/libio.h" 3 }; struct _IO_FILE { int _flags; char* _IO_read_ptr; char* _IO_read_end; char* _IO_read_base; char* _IO_write_base; char* _IO_write_ptr; char* _IO_write_end; char* _IO_buf_base; char* _IO_buf_end; char *_IO_save_base; char *_IO_backup_base; char *_IO_save_end; struct _IO_marker *_markers; struct _IO_FILE *_chain; int _fileno; int _blksize; __off_t _offset; unsigned short _cur_column; char _unused; char _shortbuf[1]; _IO_lock_t *_lock; }; typedef struct _IO_FILE _IO_FILE; struct _IO_FILE_plus; extern struct _IO_FILE_plus _IO_stdin_, _IO_stdout_, _IO_stderr_; typedef struct { __ssize_t (*read) (struct _IO_FILE *, void *, __ssize_t ) ; __ssize_t (*write) (struct _IO_FILE *, const void *, __ssize_t ) ; __off_t (*seek) (struct _IO_FILE *, __off_t , int) ; int (*close) (struct _IO_FILE *) ; } _IO_cookie_io_functions_t; struct _IO_cookie_file { struct _IO_FILE file; const void *vtable; void *cookie; _IO_cookie_io_functions_t io_functions; }; extern int __underflow (_IO_FILE*) ; extern int __uflow (_IO_FILE*) ; extern int __overflow (_IO_FILE*, int) ; extern int _IO_getc (_IO_FILE *__fp) ; extern int _IO_putc (int __c, _IO_FILE *__fp) ; extern int _IO_feof (_IO_FILE *__fp) ; extern int _IO_ferror (_IO_FILE *__fp) ; extern int _IO_peekc_locked (_IO_FILE *__fp) ; extern void _IO_flockfile (_IO_FILE *) ; extern void _IO_funlockfile (_IO_FILE *) ; extern int _IO_ftrylockfile (_IO_FILE *) ; extern int _IO_vfscanf (_IO_FILE*, const char*, __gnuc_va_list , int*) ; extern int _IO_vfprintf (_IO_FILE*, const char*, __gnuc_va_list ) ; extern __ssize_t _IO_padn (_IO_FILE *, int, __ssize_t ) ; extern size_t _IO_sgetn (_IO_FILE *, void*, size_t ) ; extern __off_t _IO_seekoff (_IO_FILE*, __off_t , int, int) ; extern __off_t _IO_seekpos (_IO_FILE*, __off_t , int) ; extern void _IO_free_backup_area (_IO_FILE*) ; # 29 "/usr/include/stdio.h" 2 3 # 69 "/usr/include/stdio.h" 3 typedef struct _IO_FILE FILE; typedef __off_t fpos_t; # 1 "/usr/include/stdio_lim.h" 1 3 # 74 "/usr/include/stdio.h" 2 3 extern FILE *stdin, *stdout, *stderr; extern void clearerr (FILE*) ; extern int fclose (FILE*) ; extern int feof (FILE*) ; extern int ferror (FILE*) ; extern int fflush (FILE*) ; extern int fgetc (FILE *) ; extern int fgetpos (FILE* fp, fpos_t *__pos) ; extern char* fgets (char*, int, FILE*) ; extern FILE* fopen (const char*, const char*) ; extern FILE* fopencookie (void *__cookie, const char *__mode, _IO_cookie_io_functions_t __io_functions) ; extern int fprintf (FILE*, const char* __format, ...) ; extern int fputc (int, FILE*) ; extern int fputs (const char *__str, FILE *__fp) ; extern size_t fread (void*, size_t, size_t, FILE*) ; extern FILE* freopen (const char*, const char*, FILE*) ; extern int fscanf (FILE *__fp, const char* __format, ...) ; extern int fseek (FILE* __fp, long int __offset, int __whence) ; extern int fsetpos (FILE* __fp, const fpos_t *__pos) ; extern long int ftell (FILE* __fp) ; extern size_t fwrite (const void*, size_t, size_t, FILE*) ; extern int getc (FILE *) ; extern int getchar (void) ; extern char* gets (char*) ; extern void perror (const char *) ; extern int printf (const char* __format, ...) ; extern int putc (int, FILE *) ; extern int putchar (int) ; extern int puts (const char *__str) ; extern int remove (const char*) ; extern int rename (const char* __old, const char* __new) ; extern void rewind (FILE*) ; extern int scanf (const char* format, ...) ; extern void setbuf (FILE*, char*) ; extern void setlinebuf (FILE*) ; extern void setbuffer (FILE*, char*, int) ; extern int setvbuf (FILE*, char*, int __mode, size_t __size) ; extern int sprintf (char*, const char* __format, ...) ; extern int sscanf (const char* string, const char* __format, ...) ; extern FILE* tmpfile (void) ; extern char* tmpnam (char*) ; extern char* tmpnam_r (char*) ; extern char *tempnam (const char *__dir, const char *__pfx) ; extern char *__stdio_gen_tempname (char *__buf, size_t bufsize, const char *dir, const char *pfx, int dir_search, size_t *lenptr, FILE **streamptr) ; extern int ungetc (int c, FILE* fp) ; extern int vfprintf (FILE *fp, char const *fmt0, __gnuc_va_list ) ; extern int vprintf (char const *fmt, __gnuc_va_list ) ; extern int vsprintf (char* string, const char* format, __gnuc_va_list ) ; extern void __libc_fatal (const char *__message) ; extern int dprintf (int, const char *, ...) ; extern int vdprintf (int, const char *, __gnuc_va_list ) ; extern int vfscanf (FILE*, const char *, __gnuc_va_list ) ; extern int __vfscanf (FILE*, const char *, __gnuc_va_list ) ; extern int vscanf (const char *, __gnuc_va_list ) ; extern int vsscanf (const char *, const char *, __gnuc_va_list ) ; extern int __vsscanf (const char *, const char *, __gnuc_va_list ) ; # 177 "/usr/include/stdio.h" 3 extern FILE *fdopen (int, const char *) ; extern int fileno (FILE*) ; extern FILE* popen (const char*, const char*) ; extern int pclose (FILE*) ; extern char *ctermid (char *__buf) ; extern char *cuserid (char * __buf) ; extern int snprintf (char *, size_t, const char *, ...) ; extern int __snprintf (char *, size_t, const char *, ...) ; extern int vsnprintf (char *, size_t, const char *, __gnuc_va_list ) ; extern int __vsnprintf (char *, size_t, const char *, __gnuc_va_list ) ; # 214 "/usr/include/stdio.h" 3 extern int __underflow (struct _IO_FILE*) ; extern int __overflow (struct _IO_FILE*, int) ; extern int sys_nerr; extern const char *const sys_errlist[]; extern void clearerr_locked (FILE *) ; extern void clearerr_unlocked (FILE *) ; extern int feof_locked (FILE *) ; extern int feof_unlocked (FILE *) ; extern int ferror_locked (FILE*) ; extern int ferror_unlocked (FILE*) ; extern int fileno_locked (FILE *) ; extern int fileno_unlocked (FILE *) ; extern int fclose_unlocked (FILE *) ; extern int fflush_locked (FILE *) ; extern int fflush_unlocked (FILE *) ; extern size_t fread_unlocked (void *, size_t, size_t, FILE *) ; extern size_t fwrite_unlocked (const void *, size_t, size_t, FILE *) ; extern int fputc_locked (int, FILE*) ; extern int fputc_unlocked (int, FILE*) ; extern int getc_locked (FILE *) ; extern int getchar_locked (void) ; extern int putc_locked (int, FILE *) ; extern int putchar_locked (int) ; extern void flockfile (FILE *) ; extern void funlockfile (FILE *) ; extern int ftrylockfile (FILE *) ; extern int getc_unlocked (FILE *) ; extern int getchar_unlocked (void) ; extern int putc_unlocked (int, FILE *) ; extern int putchar_unlocked (int) ; # 31 "f.c" 2 # 1 "/usr/include/math.h" 1 3 # 1 "/usr/include/huge_val.h" 1 3 static union { unsigned char __c[8]; double __d; } __huge_val = { { 0, 0, 0, 0, 0, 0, 0xf0, 0x7f } }; # 68 "/usr/include/huge_val.h" 3 # 33 "/usr/include/math.h" 2 3 # 1 "/usr/include/mathcalls.h" 1 3 extern double acos (double __x) ; extern double __acos (double __x) ; extern double asin (double __x) ; extern double __asin (double __x) ; extern double atan (double __x) ; extern double __atan (double __x) ; extern double atan2 (double __y, double __x) ; extern double __atan2 (double __y, double __x) ; extern double cos (double __x) ; extern double __cos (double __x) ; extern double sin (double __x) ; extern double __sin (double __x) ; extern double tan (double __x) ; extern double __tan (double __x) ; extern double cosh (double __x) ; extern double __cosh (double __x) ; extern double sinh (double __x) ; extern double __sinh (double __x) ; extern double tanh (double __x) ; extern double __tanh (double __x) ; extern double acosh (double __x) ; extern double __acosh (double __x) ; extern double asinh (double __x) ; extern double __asinh (double __x) ; extern double atanh (double __x) ; extern double __atanh (double __x) ; extern double exp (double __x) ; extern double __exp (double __x) ; extern double frexp (double __x, int *__exponent) ; extern double __frexp (double __x, int *__exponent) ; extern double ldexp (double __x, int __exponent) ; extern double __ldexp (double __x, int __exponent) ; extern double log (double __x) ; extern double __log (double __x) ; extern double log10 (double __x) ; extern double __log10 (double __x) ; extern double expm1 (double __x) ; extern double __expm1 (double __x) ; extern double log1p (double __x) ; extern double __log1p (double __x) ; extern double logb (double __x) ; extern double __logb (double __x) ; extern double modf (double __x, double *__iptr) ; extern double __modf (double __x, double *__iptr) ; extern double pow (double __x, double __y) ; extern double __pow (double __x, double __y) ; extern double sqrt (double __x) ; extern double __sqrt (double __x) ; extern double cbrt (double __x) ; extern double __cbrt (double __x) ; extern double ceil (double __x) ; extern double __ceil (double __x) ; extern double fabs (double __x) ; extern double __fabs (double __x) ; extern double floor (double __x) ; extern double __floor (double __x) ; extern double fmod (double __x, double __y) ; extern double __fmod (double __x, double __y) ; extern int isinf (double __value) ; extern int __isinf (double __value) ; extern int finite (double __value) ; extern int __finite (double __value) ; extern double copysign (double __x, double __y) ; extern double __copysign (double __x, double __y) ; extern double scalbn (double __x, int __n) ; extern double __scalbn (double __x, int __n) ; extern double drem (double __x, double __y) ; extern double __drem (double __x, double __y) ; extern double significand (double __x) ; extern double __significand (double __x) ; extern int isnan (double __value) ; extern int __isnan (double __value) ; extern int ilogb (double __x) ; extern int __ilogb (double __x) ; extern double hypot (double __x, double __y) ; extern double __hypot (double __x, double __y) ; extern double erf (double ) ; extern double __erf (double ) ; extern double erfc (double ) ; extern double __erfc (double ) ; extern double gamma (double ) ; extern double __gamma (double ) ; extern double j0 (double ) ; extern double __j0 (double ) ; extern double j1 (double ) ; extern double __j1 (double ) ; extern double jn (int, double ) ; extern double __jn (int, double ) ; extern double lgamma (double ) ; extern double __lgamma (double ) ; extern double y0 (double ) ; extern double __y0 (double ) ; extern double y1 (double ) ; extern double __y1 (double ) ; extern double yn (int, double ) ; extern double __yn (int, double ) ; extern double gamma_r (double , int *) ; extern double __gamma_r (double , int *) ; extern double lgamma_r (double , int *) ; extern double __lgamma_r (double , int *) ; extern double rint (double __x) ; extern double __rint (double __x) ; extern double scalb (double __x, double __n) ; extern double __scalb (double __x, double __n) ; extern double nextafter (double __x, double __y) ; extern double __nextafter (double __x, double __y) ; extern double remainder (double __x, double __y) ; extern double __remainder (double __x, double __y) ; # 56 "/usr/include/math.h" 2 3 # 1 "/usr/include/mathcalls.h" 1 3 extern float acosf (float __x) ; extern float __acosf (float __x) ; extern float asinf (float __x) ; extern float __asinf (float __x) ; extern float atanf (float __x) ; extern float __atanf (float __x) ; extern float atan2f (float __y, float __x) ; extern float __atan2f (float __y, float __x) ; extern float cosf (float __x) ; extern float __cosf (float __x) ; extern float sinf (float __x) ; extern float __sinf (float __x) ; extern float tanf (float __x) ; extern float __tanf (float __x) ; extern float coshf (float __x) ; extern float __coshf (float __x) ; extern float sinhf (float __x) ; extern float __sinhf (float __x) ; extern float tanhf (float __x) ; extern float __tanhf (float __x) ; extern float acoshf (float __x) ; extern float __acoshf (float __x) ; extern float asinhf (float __x) ; extern float __asinhf (float __x) ; extern float atanhf (float __x) ; extern float __atanhf (float __x) ; extern float expf (float __x) ; extern float __expf (float __x) ; extern float frexpf (float __x, int *__exponent) ; extern float __frexpf (float __x, int *__exponent) ; extern float ldexpf (float __x, int __exponent) ; extern float __ldexpf (float __x, int __exponent) ; extern float logf (float __x) ; extern float __logf (float __x) ; extern float log10f (float __x) ; extern float __log10f (float __x) ; extern float expm1f (float __x) ; extern float __expm1f (float __x) ; extern float log1pf (float __x) ; extern float __log1pf (float __x) ; extern float logbf (float __x) ; extern float __logbf (float __x) ; extern float modff (float __x, float *__iptr) ; extern float __modff (float __x, float *__iptr) ; extern float powf (float __x, float __y) ; extern float __powf (float __x, float __y) ; extern float sqrtf (float __x) ; extern float __sqrtf (float __x) ; extern float cbrtf (float __x) ; extern float __cbrtf (float __x) ; extern float ceilf (float __x) ; extern float __ceilf (float __x) ; extern float fabsf (float __x) ; extern float __fabsf (float __x) ; extern float floorf (float __x) ; extern float __floorf (float __x) ; extern float fmodf (float __x, float __y) ; extern float __fmodf (float __x, float __y) ; extern int isinff (float __value) ; extern int __isinff (float __value) ; extern int finitef (float __value) ; extern int __finitef (float __value) ; extern float copysignf (float __x, float __y) ; extern float __copysignf (float __x, float __y) ; extern float scalbnf (float __x, int __n) ; extern float __scalbnf (float __x, int __n) ; extern float dremf (float __x, float __y) ; extern float __dremf (float __x, float __y) ; extern float significandf (float __x) ; extern float __significandf (float __x) ; extern int isnanf (float __value) ; extern int __isnanf (float __value) ; extern int ilogbf (float __x) ; extern int __ilogbf (float __x) ; extern float hypotf (float __x, float __y) ; extern float __hypotf (float __x, float __y) ; extern float erff (float ) ; extern float __erff (float ) ; extern float erfcf (float ) ; extern float __erfcf (float ) ; extern float gammaf (float ) ; extern float __gammaf (float ) ; extern float j0f (float ) ; extern float __j0f (float ) ; extern float j1f (float ) ; extern float __j1f (float ) ; extern float jnf (int, float ) ; extern float __jnf (int, float ) ; extern float lgammaf (float ) ; extern float __lgammaf (float ) ; extern float y0f (float ) ; extern float __y0f (float ) ; extern float y1f (float ) ; extern float __y1f (float ) ; extern float ynf (int, float ) ; extern float __ynf (int, float ) ; extern float gammaf_r (float , int *) ; extern float __gammaf_r (float , int *) ; extern float lgammaf_r (float , int *) ; extern float __lgammaf_r (float , int *) ; extern float rintf (float __x) ; extern float __rintf (float __x) ; extern float scalbf (float __x, float __n) ; extern float __scalbf (float __x, float __n) ; extern float nextafterf (float __x, float __y) ; extern float __nextafterf (float __x, float __y) ; extern float remainderf (float __x, float __y) ; extern float __remainderf (float __x, float __y) ; # 75 "/usr/include/math.h" 2 3 # 1 "/usr/include/mathcalls.h" 1 3 extern long double acosl (long double __x) ; extern long double __acosl (long double __x) ; extern long double asinl (long double __x) ; extern long double __asinl (long double __x) ; extern long double atanl (long double __x) ; extern long double __atanl (long double __x) ; extern long double atan2l (long double __y, long double __x) ; extern long double __atan2l (long double __y, long double __x) ; extern long double cosl (long double __x) ; extern long double __cosl (long double __x) ; extern long double sinl (long double __x) ; extern long double __sinl (long double __x) ; extern long double tanl (long double __x) ; extern long double __tanl (long double __x) ; extern long double coshl (long double __x) ; extern long double __coshl (long double __x) ; extern long double sinhl (long double __x) ; extern long double __sinhl (long double __x) ; extern long double tanhl (long double __x) ; extern long double __tanhl (long double __x) ; extern long double acoshl (long double __x) ; extern long double __acoshl (long double __x) ; extern long double asinhl (long double __x) ; extern long double __asinhl (long double __x) ; extern long double atanhl (long double __x) ; extern long double __atanhl (long double __x) ; extern long double expl (long double __x) ; extern long double __expl (long double __x) ; extern long double frexpl (long double __x, int *__exponent) ; extern long double __frexpl (long double __x, int *__exponent) ; extern long double ldexpl (long double __x, int __exponent) ; extern long double __ldexpl (long double __x, int __exponent) ; extern long double logl (long double __x) ; extern long double __logl (long double __x) ; extern long double log10l (long double __x) ; extern long double __log10l (long double __x) ; extern long double expm1l (long double __x) ; extern long double __expm1l (long double __x) ; extern long double log1pl (long double __x) ; extern long double __log1pl (long double __x) ; extern long double logbl (long double __x) ; extern long double __logbl (long double __x) ; extern long double modfl (long double __x, long double *__iptr) ; extern long double __modfl (long double __x, long double *__iptr) ; extern long double powl (long double __x, long double __y) ; extern long double __powl (long double __x, long double __y) ; extern long double sqrtl (long double __x) ; extern long double __sqrtl (long double __x) ; extern long double cbrtl (long double __x) ; extern long double __cbrtl (long double __x) ; extern long double ceill (long double __x) ; extern long double __ceill (long double __x) ; extern long double fabsl (long double __x) ; extern long double __fabsl (long double __x) ; extern long double floorl (long double __x) ; extern long double __floorl (long double __x) ; extern long double fmodl (long double __x, long double __y) ; extern long double __fmodl (long double __x, long double __y) ; extern int isinfl (long double __value) ; extern int __isinfl (long double __value) ; extern int finitel (long double __value) ; extern int __finitel (long double __value) ; extern long double copysignl (long double __x, long double __y) ; extern long double __copysignl (long double __x, long double __y) ; extern long double scalbnl (long double __x, int __n) ; extern long double __scalbnl (long double __x, int __n) ; extern long double dreml (long double __x, long double __y) ; extern long double __dreml (long double __x, long double __y) ; extern long double significandl (long double __x) ; extern long double __significandl (long double __x) ; extern int isnanl (long double __value) ; extern int __isnanl (long double __value) ; extern int ilogbl (long double __x) ; extern int __ilogbl (long double __x) ; extern long double hypotl (long double __x, long double __y) ; extern long double __hypotl (long double __x, long double __y) ; extern long double erfl (long double ) ; extern long double __erfl (long double ) ; extern long double erfcl (long double ) ; extern long double __erfcl (long double ) ; extern long double gammal (long double ) ; extern long double __gammal (long double ) ; extern long double j0l (long double ) ; extern long double __j0l (long double ) ; extern long double j1l (long double ) ; extern long double __j1l (long double ) ; extern long double jnl (int, long double ) ; extern long double __jnl (int, long double ) ; extern long double lgammal (long double ) ; extern long double __lgammal (long double ) ; extern long double y0l (long double ) ; extern long double __y0l (long double ) ; extern long double y1l (long double ) ; extern long double __y1l (long double ) ; extern long double ynl (int, long double ) ; extern long double __ynl (int, long double ) ; extern long double gammal_r (long double , int *) ; extern long double __gammal_r (long double , int *) ; extern long double lgammal_r (long double , int *) ; extern long double __lgammal_r (long double , int *) ; extern long double rintl (long double __x) ; extern long double __rintl (long double __x) ; extern long double scalbl (long double __x, long double __n) ; extern long double __scalbl (long double __x, long double __n) ; extern long double nextafterl (long double __x, long double __y) ; extern long double __nextafterl (long double __x, long double __y) ; extern long double remainderl (long double __x, long double __y) ; extern long double __remainderl (long double __x, long double __y) ; # 91 "/usr/include/math.h" 2 3 extern int signgam; typedef enum { _IEEE_ = -1, _SVID_, _XOPEN_, _POSIX_ } _LIB_VERSION_TYPE; extern _LIB_VERSION_TYPE _LIB_VERSION; struct exception { int type; char *name; double arg1; double arg2; double retval; }; extern int __matherr (struct exception *) ; extern int matherr (struct exception *) ; # 1 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/float.h" 1 3 union __convert_long_double { unsigned __convert_long_double_i[4]; long double __convert_long_double_d; }; # 155 "/usr/include/math.h" 2 3 # 165 "/usr/include/math.h" 3 # 202 "/usr/include/math.h" 3 # 32 "f.c" 2 # 1 "/usr/include/string.h" 1 3 # 1 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 1 3 # 19 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 61 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 131 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 182 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 258 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 270 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 302 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 33 "/usr/include/string.h" 2 3 extern void * memcpy (void * __dest, const void * __src, size_t __n) ; extern void * memmove (void * __dest, const void * __src, size_t __n) ; extern void * __memccpy (void * __dest, const void * __src, int __c, size_t __n) ; extern void * memccpy (void * __dest, const void * __src, int __c, size_t __n) ; extern void * memset (void * __s, int __c, size_t __n) ; extern int memcmp (const void * __s1, const void * __s2, size_t __n) ; extern void * memchr (const void * __s, int __c, size_t __n) ; extern char *strcpy (char *__dest, const char *__src) ; extern char *strncpy (char *__dest, const char *__src, size_t __n) ; extern char *strcat (char *__dest, const char *__src) ; extern char *strncat (char *__dest, const char *__src, size_t __n) ; extern int strcmp (const char *__s1, const char *__s2) ; extern int strncmp (const char *__s1, const char *__s2, size_t __n) ; extern int strcoll (const char *__s1, const char *__s2) ; extern size_t strxfrm (char *__dest, const char *__src, size_t __n) ; extern char *__strdup (const char *__s) ; extern char *strdup (const char *__s) ; extern char *__strndup (const char *__string, size_t __n) ; # 121 "/usr/include/string.h" 3 extern char *strchr (const char *__s, int __c) ; extern char *strrchr (const char *__s, int __c) ; extern size_t strcspn (const char *__s, const char *__reject) ; extern size_t strspn (const char *__s, const char *__accept) ; extern char *strpbrk (const char *__s, const char *__accept) ; extern char *strstr (const char *__haystack, const char *__needle) ; extern char *strtok (char *__s, const char *__delim) ; extern char *strtok_r (char *__s, const char *__delim, char **__save_ptr) ; extern size_t strlen (const char *__s) ; # 173 "/usr/include/string.h" 3 extern char *strerror (int __errnum) ; extern char *__strerror_r (int __errnum, char *__buf, size_t __buflen) ; extern char *strerror_r (int __errnum, char *__buf, size_t __buflen) ; extern void bcopy (const void * __src, void * __dest, size_t __n) ; extern void bzero (void * __s, size_t __n) ; extern int bcmp (const void * __s1, const void * __s2, size_t __n) ; extern char *index (const char *__s, int __c) ; extern char *rindex (const char *__s, int __c) ; extern int ffs (int __i) ; extern int __strcasecmp (const char *__s1, const char *__s2) ; extern int strcasecmp (const char *__s1, const char *__s2) ; extern int __strncasecmp (const char *__s1, const char *__s2, size_t __n) ; extern int strncasecmp (const char *__s1, const char *__s2, size_t __n) ; extern char *__strsep (char **__stringp, const char *__delim) ; extern char *strsep (char **__stringp, const char *__delim) ; # 241 "/usr/include/string.h" 3 extern char *basename (const char *__filename) ; # 33 "f.c" 2 # 1 "/usr/include/unistd.h" 1 3 # 1 "/usr/include/posix_opt.h" 1 3 # 137 "/usr/include/unistd.h" 2 3 typedef __ssize_t ssize_t; # 1 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 1 3 # 19 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 61 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 131 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 182 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 258 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 270 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 302 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 157 "/usr/include/unistd.h" 2 3 extern int __access (const char *__name, int __type) ; extern int access (const char *__name, int __type) ; extern __off_t __lseek (int __fd, __off_t __offset, int __whence) ; extern __off_t lseek (int __fd, __off_t __offset, int __whence) ; extern int __close (int __fd) ; extern int close (int __fd) ; extern ssize_t __read (int __fd, void * __buf, size_t __nbytes) ; extern ssize_t read (int __fd, void * __buf, size_t __nbytes) ; extern ssize_t __write (int __fd, const void * __buf, size_t __n) ; extern ssize_t write (int __fd, const void * __buf, size_t __n) ; extern int __pipe (int __pipedes[2]) ; extern int pipe (int __pipedes[2]) ; extern unsigned int alarm (unsigned int __seconds) ; extern unsigned int sleep (unsigned int __seconds) ; extern unsigned int ualarm (unsigned int __value, unsigned int __interval) ; extern void usleep (unsigned int __useconds) ; extern int pause (void) ; extern int __chown (const char *__file, __uid_t __owner, __gid_t __group) ; extern int chown (const char *__file, __uid_t __owner, __gid_t __group) ; extern int __fchown (int __fd, __uid_t __owner, __gid_t __group) ; extern int fchown (int __fd, __uid_t __owner, __gid_t __group) ; extern int __lchown (const char *__file, __uid_t __owner, __gid_t __group) ; extern int lchown (const char *__file, __uid_t __owner, __gid_t __group) ; extern int __chdir (const char *__path) ; extern int chdir (const char *__path) ; extern int fchdir (int __fd) ; extern char *__getcwd (char *__buf, size_t __size) ; extern char *getcwd (char *__buf, size_t __size) ; # 317 "/usr/include/unistd.h" 3 extern char *getwd (char *__buf) ; extern int __dup (int __fd) ; extern int dup (int __fd) ; extern int __dup2 (int __fd, int __fd2) ; extern int dup2 (int __fd, int __fd2) ; extern char **__environ; extern int __execve (const char *__path, char * const __argv[], char * const __envp[]) ; extern int execve (const char *__path, char * const __argv[], char * const __envp[]) ; extern int execv (const char *__path, char * const __argv[]) ; extern int execle (const char *__path, const char *__arg, ...) ; extern int execl (const char *__path, const char *__arg, ...) ; extern int execvp (const char *__file, char * const __argv[]) ; extern int execlp (const char *__file, const char *__arg, ...) ; extern int nice (int __inc) ; extern void _exit (int __status) ; # 1 "/usr/include/confname.h" 1 3 enum { _PC_LINK_MAX, _PC_MAX_CANON, _PC_MAX_INPUT, _PC_NAME_MAX, _PC_PATH_MAX, _PC_PIPE_BUF, _PC_CHOWN_RESTRICTED, _PC_NO_TRUNC, _PC_VDISABLE, _PC_SYNC_IO, _PC_ASYNC_IO, _PC_PRIO_IO, _PC_SOCK_MAXBUF }; enum { _SC_ARG_MAX, _SC_CHILD_MAX, _SC_CLK_TCK, _SC_NGROUPS_MAX, _SC_OPEN_MAX, _SC_STREAM_MAX, _SC_TZNAME_MAX, _SC_JOB_CONTROL, _SC_SAVED_IDS, _SC_REALTIME_SIGNALS, _SC_PRIORITY_SCHEDULING, _SC_TIMERS, _SC_ASYNCHRONOUS_IO, _SC_PRIORITIZED_IO, _SC_SYNCHRONIZED_IO, _SC_FSYNC, _SC_MAPPED_FILES, _SC_MEMLOCK, _SC_MEMLOCK_RANGE, _SC_MEMORY_PROTECTION, _SC_MESSAGE_PASSING, _SC_SEMAPHORES, _SC_SHARED_MEMORY_OBJECTS, _SC_AIO_LISTIO_MAX, _SC_AIO_MAX, _SC_AIO_PRIO_DELTA_MAX, _SC_DELAYTIMER_MAX, _SC_MQ_OPEN_MAX, _SC_MQ_PRIO_MAX, _SC_VERSION, _SC_PAGESIZE, _SC_RTSIG_MAX, _SC_SEM_NSEMS_MAX, _SC_SEM_VALUE_MAX, _SC_SIGQUEUE_MAX, _SC_TIMER_MAX, _SC_BC_BASE_MAX, _SC_BC_DIM_MAX, _SC_BC_SCALE_MAX, _SC_BC_STRING_MAX, _SC_COLL_WEIGHTS_MAX, _SC_EQUIV_CLASS_MAX, _SC_EXPR_NEST_MAX, _SC_LINE_MAX, _SC_RE_DUP_MAX, _SC_CHARCLASS_NAME_MAX, _SC_2_VERSION, _SC_2_C_BIND, _SC_2_C_DEV, _SC_2_FORT_DEV, _SC_2_FORT_RUN, _SC_2_SW_DEV, _SC_2_LOCALEDEF, _SC_PII, _SC_PII_XTI, _SC_PII_SOCKET, _SC_PII_INTERNET, _SC_PII_OSI, _SC_POLL, _SC_SELECT, _SC_UIO_MAXIOV, _SC_PII_INTERNET_STREAM, _SC_PII_INTERNET_DGRAM, _SC_PII_OSI_COTS, _SC_PII_OSI_CLTS, _SC_PII_OSI_M, _SC_T_IOV_MAX, _SC_THREADS, _SC_THREAD_SAFE_FUNCTIONS, _SC_GETGR_R_SIZE_MAX, _SC_GETPW_R_SIZE_MAX, _SC_LOGIN_NAME_MAX, _SC_TTY_NAME_MAX, _SC_THREAD_DESTRUCTOR_ITERATIONS, _SC_THREAD_KEYS_MAX, _SC_THREAD_STACK_MIN, _SC_THREAD_THREADS_MAX, _SC_THREAD_ATTR_STACKADDR, _SC_THREAD_ATTR_STACKSIZE, _SC_THREAD_PRIORITY_SCHEDULING, _SC_THREAD_PRIO_INHERIT, _SC_THREAD_PRIO_PROTECT, _SC_THREAD_PROCESS_SHARED, _SC_NPROCESSORS_CONF, _SC_NPROCESSORS_ONLN, _SC_PHYS_PAGES, _SC_AVPHYS_PAGES, _SC_ATEXIT_MAX, _SC_PASS_MAX, _SC_XOPEN_VERSION, _SC_XOPEN_XCU_VERSION, _SC_XOPEN_UNIX, _SC_XOPEN_CRYPT, _SC_XOPEN_ENH_I18N, _SC_XOPEN_SHM, _SC_2_CHAR_TERM, _SC_2_C_VERSION, _SC_2_UPE, _SC_XOPEN_XPG2, _SC_XOPEN_XPG3, _SC_XOPEN_XPG4, _SC_CHAR_BIT, _SC_CHAR_MAX, _SC_CHAR_MIN, _SC_INT_MAX, _SC_INT_MIN, _SC_LONG_BIT, _SC_WORD_BIT, _SC_MB_LEN_MAX, _SC_NZERO, _SC_SSIZE_MAX, _SC_SCHAR_MAX, _SC_SCHAR_MIN, _SC_SHRT_MAX, _SC_SHRT_MIN, _SC_UCHAR_MAX, _SC_UINT_MAX, _SC_ULONG_MAX, _SC_USHRT_MAX, _SC_NL_ARGMAX, _SC_NL_LANGMAX, _SC_NL_MSGMAX, _SC_NL_NMAX, _SC_NL_SETMAX, _SC_NL_TEXTMAX }; enum { _CS_PATH }; # 392 "/usr/include/unistd.h" 2 3 extern long int __pathconf (const char *__path, int __name) ; extern long int pathconf (const char *__path, int __name) ; extern long int __fpathconf (int __fd, int __name) ; extern long int fpathconf (int __fd, int __name) ; extern long int __sysconf (int __name) ; extern long int sysconf (int __name) ; extern size_t confstr (int __name, char *__buf, size_t __len) ; extern __pid_t __getpid (void) ; extern __pid_t getpid (void) ; extern __pid_t __getppid (void) ; extern __pid_t getppid (void) ; extern __pid_t getpgrp (void) ; extern int setpgid (__pid_t __pid, __pid_t __pgid) ; extern __pid_t __getpgid (__pid_t __pid) ; extern int setpgrp (void) ; extern __pid_t __setsid (void) ; extern __pid_t setsid (void) ; extern __uid_t __getuid (void) ; extern __uid_t getuid (void) ; extern __uid_t __geteuid (void) ; extern __uid_t geteuid (void) ; extern __gid_t __getgid (void) ; extern __gid_t getgid (void) ; extern __gid_t __getegid (void) ; extern __gid_t getegid (void) ; extern int __getgroups (int __size, __gid_t __list[]) ; extern int getgroups (int __size, __gid_t __list[]) ; extern int __setuid (__uid_t __uid) ; extern int setuid (__uid_t __uid) ; extern int __setreuid (__uid_t __ruid, __uid_t __euid) ; extern int setreuid (__uid_t __ruid, __uid_t __euid) ; extern int seteuid (__uid_t __uid) ; extern int __setgid (__gid_t __gid) ; extern int setgid (__gid_t __gid) ; extern int __setregid (__gid_t __rgid, __gid_t __egid) ; extern int setregid (__gid_t __rgid, __gid_t __egid) ; extern int setegid (__gid_t __gid) ; extern __pid_t __fork (void) ; extern __pid_t fork (void) ; extern __pid_t __vfork (void) ; extern __pid_t vfork (void) ; extern char *ttyname (int __fd) ; extern int __ttyname_r (int __fd, char *__buf, size_t __buflen) ; extern int ttyname_r (int __fd, char *__buf, size_t __buflen) ; extern int __isatty (int __fd) ; extern int isatty (int __fd) ; extern int ttyslot (void) ; extern int __link (const char *__from, const char *__to) ; extern int link (const char *__from, const char *__to) ; extern int __symlink (const char *__from, const char *__to) ; extern int symlink (const char *__from, const char *__to) ; extern int __readlink (const char *__path, char *__buf, size_t __len) ; extern int readlink (const char *__path, char *__buf, size_t __len) ; extern int __unlink (const char *__name) ; extern int unlink (const char *__name) ; extern int __rmdir (const char *__path) ; extern int rmdir (const char *__path) ; extern __pid_t tcgetpgrp (int __fd) ; extern int tcsetpgrp (int __fd, __pid_t __pgrp_id) ; extern char *getlogin (void) ; extern int setlogin (const char *__name) ; extern int getopt (int __argc, char * const * __argv, const char *__opts) ; extern int opterr; extern int optind; extern int optopt; extern char *optarg; extern int __gethostname (char *__name, size_t __len) ; extern int gethostname (char *__name, size_t __len) ; extern int sethostname (const char *__name, size_t __len) ; extern int sethostid (long int __id) ; extern int getdomainname (char *__name, size_t __len) ; extern int setdomainname (const char *__name, size_t __len) ; extern int fsync (int __fd) ; extern int vhangup (void) ; extern int revoke (const char *__file) ; extern int profil (unsigned short int *__sample_buffer, size_t __size, size_t __offset, unsigned int __scale) ; extern int acct (const char *__name) ; extern int chroot (const char *__path) ; extern char *getusershell (void) ; extern void endusershell (void) ; extern void setusershell (void) ; extern char *getpass (const char *__prompt) ; extern int daemon (int __nochdir, int __noclose) ; extern long int gethostid (void) ; extern int sync (void) ; extern int __getpagesize (void) ; extern int getpagesize (void) ; extern int truncate (const char *__file, __off_t __length) ; extern int ftruncate (int __fd, __off_t __length) ; extern int __getdtablesize (void) ; extern int getdtablesize (void) ; extern int __brk (void * __addr) ; extern int brk (void * __addr) ; # 1 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 1 3 # 19 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 61 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 typedef int ptrdiff_t; # 184 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 258 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 270 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 302 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 768 "/usr/include/unistd.h" 2 3 extern void * __sbrk (ptrdiff_t __delta) ; extern void * sbrk (ptrdiff_t __delta) ; extern long int syscall (long int __sysno, ...) ; extern int lockf (int __fd, int __cmd, __off_t __len) ; # 833 "/usr/include/unistd.h" 3 extern int fdatasync (int __fildes) ; # 861 "/usr/include/unistd.h" 3 # 34 "f.c" 2 # 1 "utils.h" 1 # 1 "/usr/include/time.h" 1 3 # 1 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 1 3 # 19 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 61 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 131 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 182 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 258 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 270 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 302 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 38 "/usr/include/time.h" 2 3 # 1 "/usr/include/timebits.h" 1 3 # 32 "/usr/include/timebits.h" 3 # 47 "/usr/include/time.h" 2 3 typedef __clock_t clock_t; typedef __time_t time_t; struct timespec { long int tv_sec; long int tv_nsec; }; struct tm { int tm_sec; int tm_min; int tm_hour; int tm_mday; int tm_mon; int tm_year; int tm_wday; int tm_yday; int tm_isdst; long int tm_gmtoff; const char *tm_zone; }; extern clock_t clock (void) ; extern time_t time (time_t *__timer) ; extern double difftime (time_t __time1, time_t __time0) ; extern time_t mktime (struct tm *__tp) ; extern time_t __mktime_internal (struct tm *__tp, struct tm *(*__func) (const time_t *, struct tm *), time_t *__offset) ; extern size_t strftime (char *__s, size_t __maxsize, const char *__format, const struct tm *__tp) ; extern struct tm *gmtime (const time_t *__timer) ; extern struct tm *localtime (const time_t *__timer) ; extern struct tm *__gmtime_r (const time_t *__timer, struct tm *__tp) ; extern struct tm *gmtime_r (const time_t *__timer, struct tm *__tp) ; extern struct tm *__localtime_r (const time_t *__timer, struct tm *__tp) ; extern struct tm *localtime_r (const time_t *__timer, struct tm *__tp) ; extern int __offtime (const time_t *__timer, long int __offset, struct tm *__TP) ; extern char *asctime (const struct tm *__tp) ; extern char *ctime (const time_t *__timer) ; extern char *__asctime_r (const struct tm *__tp, char *__buf) ; extern char *asctime_r (const struct tm *__tp, char *__buf) ; extern char *ctime_r (const time_t *__timer, char *__buf) ; extern char *__tzname[2]; extern int __daylight; extern long int __timezone; extern void __tzset (void) ; extern char *tzname[2]; extern long int __tzname_max (void) ; extern void tzset (void) ; extern int daylight; extern long int timezone; extern int stime (const time_t *__when) ; extern time_t timegm (struct tm *__tp) ; extern time_t timelocal (struct tm *__tp) ; extern int dysize (int __year) ; extern int __nanosleep (const struct timespec *__requested_time, struct timespec *__remaining) ; extern int nanosleep (const struct timespec *__requested_time, struct timespec *__remaining) ; # 7 "utils.h" 2 # 1 "/usr/include/stdlib.h" 1 3 # 1 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 1 3 # 19 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 61 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 131 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 182 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 typedef long int wchar_t; # 270 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 302 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 32 "/usr/include/stdlib.h" 2 3 typedef struct { int quot; int rem; } div_t; typedef struct { long int quot; long int rem; } ldiv_t; extern int __ctype_get_mb_cur_max (void) ; extern double atof (const char *__nptr) ; extern int atoi (const char *__nptr) ; extern long int atol (const char *__nptr) ; extern double strtod (const char *__nptr, char **__endptr) ; extern long int strtol (const char *__nptr, char **__endptr, int __base) ; extern unsigned long int strtoul (const char *__nptr, char **__endptr, int __base) ; # 121 "/usr/include/stdlib.h" 3 extern double __strtod_internal (const char *__nptr, char **__endptr, int __group) ; extern float __strtof_internal (const char *__nptr, char **__endptr, int __group) ; extern long double __strtold_internal (const char *__nptr, char **__endptr, int __group) ; extern long int __strtol_internal (const char *__nptr, char **__endptr, int __base, int __group) ; extern unsigned long int __strtoul_internal (const char *__nptr, char **__endptr, int __base, int __group) ; # 147 "/usr/include/stdlib.h" 3 # 197 "/usr/include/stdlib.h" 3 extern char *l64a (long int __n) ; extern long int a64l (const char *__s) ; # 1 "/usr/include/sys/types.h" 1 3 typedef __u_char u_char; typedef __u_short u_short; typedef __u_int u_int; typedef __u_long u_long; typedef __quad_t quad_t; typedef __u_quad_t u_quad_t; typedef __fsid_t fsid_t; typedef __dev_t dev_t; typedef __gid_t gid_t; typedef __ino_t ino_t; typedef __mode_t mode_t; typedef __nlink_t nlink_t; typedef __off_t off_t; typedef __loff_t loff_t; typedef __pid_t pid_t; typedef __uid_t uid_t; typedef __daddr_t daddr_t; typedef __caddr_t caddr_t; typedef __key_t key_t; # 1 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 1 3 # 19 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 61 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 131 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 182 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 258 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 270 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 302 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 70 "/usr/include/sys/types.h" 2 3 typedef unsigned long int ulong; typedef unsigned short int ushort; typedef unsigned int uint; typedef char int8_t; typedef unsigned char u_int8_t; typedef short int int16_t; typedef unsigned short int u_int16_t; typedef int int32_t; typedef unsigned int u_int32_t; typedef int register_t; # 118 "/usr/include/sys/types.h" 3 # 1 "/usr/include/endian.h" 1 3 # 1 "/usr/include/bytesex.h" 1 3 # 34 "/usr/include/endian.h" 2 3 # 123 "/usr/include/sys/types.h" 2 3 # 1 "/usr/include/sys/select.h" 1 3 # 1 "/usr/include/selectbits.h" 1 3 # 51 "/usr/include/selectbits.h" 3 # 31 "/usr/include/sys/select.h" 2 3 struct timeval; typedef __fd_mask fd_mask; typedef __fd_set fd_set; extern int __select (int __nfds, __fd_set *__readfds, __fd_set *__writefds, __fd_set *__exceptfds, struct timeval *__timeout) ; extern int select (int __nfds, __fd_set *__readfds, __fd_set *__writefds, __fd_set *__exceptfds, struct timeval *__timeout) ; extern int __pselect (int __nfds, __fd_set *__readfds, __fd_set *__writefds, __fd_set *__exceptfds, struct timespec *__timeout) ; extern int pselect (int __nfds, __fd_set *__readfds, __fd_set *__writefds, __fd_set *__exceptfds, struct timespec *__timeout) ; # 126 "/usr/include/sys/types.h" 2 3 # 210 "/usr/include/stdlib.h" 2 3 extern int32_t __random (void) ; extern int32_t random (void) ; extern void __srandom (unsigned int __seed) ; extern void srandom (unsigned int __seed) ; extern void * __initstate (unsigned int __seed, void * __statebuf, size_t __statelen) ; extern void * initstate (unsigned int __seed, void * __statebuf, size_t __statelen) ; extern void * __setstate (void * __statebuf) ; extern void * setstate (void * __statebuf) ; struct random_data { int32_t *fptr; int32_t *rptr; int32_t *state; int rand_type; int rand_deg; int rand_sep; int32_t *end_ptr; }; extern int __random_r (struct random_data *__buf, int32_t *__result) ; extern int random_r (struct random_data *__buf, int32_t *__result) ; extern int __srandom_r (unsigned int __seed, struct random_data *__buf) ; extern int srandom_r (unsigned int __seed, struct random_data *__buf) ; extern int __initstate_r (unsigned int __seed, void * __statebuf, size_t __statelen, struct random_data *__buf) ; extern int initstate_r (unsigned int __seed, void * __statebuf, size_t __statelen, struct random_data *__buf) ; extern int __setstate_r (void * __statebuf, struct random_data *__buf) ; extern int setstate_r (void * __statebuf, struct random_data *__buf) ; extern int rand (void) ; extern void srand (unsigned int __seed) ; extern int __rand_r (unsigned int *__seed) ; extern int rand_r (unsigned int *__seed) ; extern double drand48 (void) ; extern double erand48 (unsigned short int __xsubi[3]) ; extern long lrand48 (void) ; extern long nrand48 (unsigned short int __xsubi[3]) ; extern long mrand48 (void) ; extern long jrand48 (unsigned short int __xsubi[3]) ; extern void srand48 (long __seedval) ; extern unsigned short int *seed48 (unsigned short int __seed16v[3]) ; extern void lcong48 (unsigned short int __param[7]) ; struct drand48_data { unsigned short int x[3]; unsigned short int a[3]; unsigned short int c; unsigned short int old_x[3]; int init; }; extern int drand48_r (struct drand48_data *__buffer, double *__result) ; extern int erand48_r (unsigned short int __xsubi[3], struct drand48_data *__buffer, double *__result) ; extern int lrand48_r (struct drand48_data *__buffer, long *__result) ; extern int nrand48_r (unsigned short int __xsubi[3], struct drand48_data *__buffer, long *__result) ; extern int mrand48_r (struct drand48_data *__buffer, long *__result) ; extern int jrand48_r (unsigned short int __xsubi[3], struct drand48_data *__buffer, long *__result) ; extern int srand48_r (long __seedval, struct drand48_data *__buffer) ; extern int seed48_r (unsigned short int __seed16v[3], struct drand48_data *__buffer) ; extern int lcong48_r (unsigned short int __param[7], struct drand48_data *__buffer) ; extern int __drand48_iterate (unsigned short int __xsubi[3], struct drand48_data *__buffer) ; extern void * malloc (size_t __size) ; extern void * realloc (void * __ptr, size_t __size) ; extern void * calloc (size_t __nmemb, size_t __size) ; extern void free (void * __ptr) ; extern void cfree (void * __ptr) ; # 1 "/usr/include/alloca.h" 1 3 # 1 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 1 3 # 19 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 61 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 131 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 182 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 258 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 270 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 302 "/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.90.29/include/stddef.h" 3 # 25 "/usr/include/alloca.h" 2 3 extern void * __alloca (size_t __size) ; extern void * alloca (size_t __size) ; # 360 "/usr/include/stdlib.h" 2 3 extern void * valloc (size_t __size) ; extern void abort (void) ; extern int atexit (void (*__func) (void)) ; extern int __on_exit (void (*__func) (int __status, void * __arg), void * __arg) ; extern int on_exit (void (*__func) (int __status, void * __arg), void * __arg) ; extern void exit (int __status) ; extern char *getenv (const char *__name) ; extern char *__secure_getenv (const char *__name) ; extern int putenv (const char *__string) ; extern int setenv (const char *__name, const char *__value, int __replace) ; extern void unsetenv (const char *__name) ; extern int __clearenv (void) ; extern int clearenv (void) ; extern char *mktemp (char *__template) ; extern int mkstemp (char *__template) ; extern int system (const char *__command) ; extern char *realpath (const char *__name, char *__resolved) ; typedef int (*__compar_fn_t) (const void * , const void * ) ; extern void * bsearch (const void * __key, const void * __base, size_t __nmemb, size_t __size, __compar_fn_t __compar) ; extern void qsort (void * __base, size_t __nmemb, size_t __size, __compar_fn_t __compar) ; extern int abs (int __x) ; extern long int labs (long int __x) ; extern div_t div (int __numer, int __denom) ; extern ldiv_t ldiv (long int __numer, long int __denom) ; extern char *ecvt (double __value, int __ndigit, int *__decpt, int *__sign) ; extern char *fcvt (double __value, int __ndigit, int *__decpt, int *__sign) ; extern char *gcvt (double __value, int __ndigit, char *__buf) ; extern char *qecvt (long double __value, int __ndigit, int *__decpt, int *__sign) ; extern char *qfcvt (long double __value, int __ndigit, int *__decpt, int *__sign) ; extern char *qgcvt (long double __value, int __ndigit, char *__buf) ; extern int ecvt_r (double __value, int __ndigit, int *__decpt, int *__sign, char *__buf, size_t __len) ; extern int fcvt_r (double __value, int __ndigit, int *__decpt, int *__sign, char *__buf, size_t __len) ; extern int qecvt_r (long double __value, int __ndigit, int *__decpt, int *__sign, char *__buf, size_t __len) ; extern int qfcvt_r (long double __value, int __ndigit, int *__decpt, int *__sign, char *__buf, size_t __len) ; extern int mblen (const char *__s, size_t __n) ; extern int mbtowc (wchar_t *__pwc, const char *__s, size_t __n) ; extern int wctomb (char *__s, wchar_t __wchar) ; extern size_t mbstowcs (wchar_t *__pwcs, const char *__s, size_t __n) ; extern size_t wcstombs (char *__s, const wchar_t *__pwcs, size_t __n) ; extern int rpmatch (const char *__response) ; # 582 "/usr/include/stdlib.h" 3 # 9 "utils.h" 2 # 1 "/usr/include/sys/time.h" 1 3 # 1 "/usr/include/timebits.h" 1 3 struct timeval { time_t tv_sec; time_t tv_usec; }; # 50 "/usr/include/timebits.h" 3 # 28 "/usr/include/sys/time.h" 2 3 struct timezone { int tz_minuteswest; int tz_dsttime; }; extern int __gettimeofday (struct timeval *__tv, struct timezone *__tz) ; extern int gettimeofday (struct timeval *__tv, struct timezone *__tz) ; extern int __settimeofday (const struct timeval *__tv, const struct timezone *__tz) ; extern int settimeofday (const struct timeval *__tv, const struct timezone *__tz) ; extern int __adjtime (const struct timeval *__delta, struct timeval *__olddelta) ; extern int adjtime (const struct timeval *__delta, struct timeval *__olddelta) ; enum __itimer_which { ITIMER_REAL = 0, ITIMER_VIRTUAL = 1, ITIMER_PROF = 2 }; struct itimerval { struct timeval it_interval; struct timeval it_value; }; extern int __getitimer (enum __itimer_which __which, struct itimerval *__value) ; extern int getitimer (enum __itimer_which __which, struct itimerval *__value) ; extern int __setitimer (enum __itimer_which __which, const struct itimerval *__new, struct itimerval *__old) ; extern int setitimer (enum __itimer_which __which, const struct itimerval *__new, struct itimerval *__old) ; extern int __utimes (const char *__file, struct timeval __tvp[2]) ; extern int utimes (const char *__file, struct timeval __tvp[2]) ; # 145 "/usr/include/sys/time.h" 3 # 154 "/usr/include/sys/time.h" 3 # 11 "utils.h" 2 struct schedule_list { char token[32 ]; int n; double *s; int nontrivial; struct schedule_list *next; }; double frandom(); void randnorm (int n, double x[]); int freetowrite (char filename[]); double watch(char token[], int command); void waitfor (double seconds); char *earth_time (double seconds); void mateqv(double A[3][3], double B[3][3]); void matran(double A[3][3], double B[3][3]); void matadd(double A[3][3], double B[3][3], double C[3][3]); void matsub(double A[3][3], double B[3][3], double C[3][3]); void matmul(double A[3][3], double B[3][3], double C[3][3]); double matinv(double A[3][3], double B[3][3]); char *diag3 (double A[3][3], double eigval[3], double Q[3][3]); double trace_decompose (double A[3][3], double B[3][3], double C[3][3]); void quicksort (int n, double arr[], int idx[]); struct schedule_list *add_schedule (char token[], char buf[], double tmin, double tmax); double schedule (char token[], double t); double *calm_period (char token[], double t1, double t2); void cancel_schedule (char token[]); void cancel_all_schedules(); void make_space (int idx[], int list[], int i, int min, int max); void safe_append (int idx[], int list[], int value, int i, int min, int max); # 35 "f.c" 2 # 1 "smile.h" 1 static const char *fn_freq[] = {"freq.out", "shear_freq.out"}; static const char *fn_corr[] = {"corr.out", "shear_corr.out"}; static const char *fn_flux[][3] = {{"jx.out", "jy.out", "jz.out"}, {"syz.out", "sxz.out", "sxy.out"}}; double smile (int which, FILE *output); # 37 "f.c" 2 static int LJ6_12, WOON, CALC, NTH, NTS, HSCH, SSCH, TSCH, VOLUMETRIC, HCALM, SCALM, TCALM, USE_CONFIG_H, maxkb, kstart, kw_lp, iseed, np, np3, neff, kw_gr, nw_repartition, nw_anchor, anchortum, kb, nkind, ntpd, nk, nar, scale_H, scale_T, Havg_start, nmb, *mb, *tpd, *mybin, *binbindx, *binbinlist, *bindx, *binlist, *idx, *list, nc[4], maxbin[4], kind[128 ][2]; static double MAX_STRAIN, MAX_DRIFT, RLIST, GRCUT, TR, delta, tau, wallmass, volume, volumeold, hydro, shear, T, dTdR, totalmass, r2del, kinewall, HTnow, msd, msd_kstart, structure_factor, pote, kine, Tnow, dE2, pressure, hamiltonian, factorT, factorHT, Treal, dV2, dEdV, bulk_modulus, alpha, Cv, Cp, *sa, *s, *s1, *s2, *s3, *s4, *s5, *v, *f, *fs, *d, *mass, *ep, l_ratio[3], e_ratio[3], maxl_ratio, minl_ratio[2], r[4], x[12], cj[3], c[6][6], H[3][3], H1[3][3], H2[3][3], H3[3][3], H4[3][3], H5[3][3], HI[3][3], HIOLD[3][3], A[3][3], B[3][3], G[3][3], GI[3][3], GA[3][3], G1[3][3], HT[3][3], H1T[3][3], HIT[3][3], Havg[3][3], sout[3][3], virial[3][3], stress[3][3], fH[3][3], HIBEGIN[3][3]; static int nvoids, tmp_config_freq; static double voids[128 ][4], df[4], af[5], density[100 ], vx[100 ], temp[100 ], taum[100 ][3][3]; static char tmp_config[128 ]; static double kinetum, k1tum, k2tum, Htum[3][3], strestum[3][3], v1k1tum[3][3], v1k2tum[3][3], v2k2tum[3][3][3][3], ck1tum[6][6]; static double kinesum, k1sum, k2sum, potesum, Hsum[3][3], stressum[3][3], v1k1sum[3][3], v1k2sum[3][3], v2k2sum[3][3][3][3], ck1sum[6][6]; static char buf[512], fn_config_read[128 ], **tp, fn_config_write[128 ], potfun_name[128 ], fn_property[] = "temp.out", fn_structure[] = "struc.out", *Hschname[6] = {"H11","H12","H13","H22","H23","H33"}, *Sschname[6] = {"S11","S12","S13","S22","S23","S33"}; static FILE *config,*property,*structure,*heat_current[3],*shear_stress[3]; static void write_config (char filename[]); static void save_tmp_config (); static void save_mesh (char filename[]); static void ds (double s1[], double s2[], double d[]); static double r2 (int i, int j, double si[], double sj[]); static void pair_potential(); static double T_real_to_T_md (double TR, double *dTdR); static double T_md_to_T_real (double T); static void divide_supercell_into_bins(); static void particles_move_in(); static void transfer_particle (int i, int o, int n); static void re_anchor (int i); static void initialize_all_lists(); static void print_list_statistics (FILE *out); static void accumulate_and_save_gr (); static void print_atom_statistics (FILE *out, char *sp); static int TFE (double sz1, double sz2); void write_config (char filename[]) { int i; FILE *config; config = fopen (filename, "w+"); fprintf (config, "Number of particles = %d\n\n", np); fprintf (config, "H(1,1) = %f A\n", H[0][0]* ((3.405E-10) *1E10) ); fprintf (config, "H(1,2) = %f A\n", H[0][1]* ((3.405E-10) *1E10) ); fprintf (config, "H(1,3) = %f A\n", H[0][2]* ((3.405E-10) *1E10) ); fprintf (config, "nc(1) = %d\n\n", nc[0]); fprintf (config, "H(2,1) = %f A\n", H[1][0]* ((3.405E-10) *1E10) ); fprintf (config, "H(2,2) = %f A\n", H[1][1]* ((3.405E-10) *1E10) ); fprintf (config, "H(2,3) = %f A\n", H[1][2]* ((3.405E-10) *1E10) ); fprintf (config, "nc(2) = %d\n\n", nc[1]); fprintf (config, "H(3,1) = %f A\n", H[2][0]* ((3.405E-10) *1E10) ); fprintf (config, "H(3,2) = %f A\n", H[2][1]* ((3.405E-10) *1E10) ); fprintf (config, "H(3,3) = %f A\n", H[2][2]* ((3.405E-10) *1E10) ); fprintf (config, "nc(3) = %d\n\n", nc[2]); fprintf (config, "TP Mass(amu) sx sy sz "); fprintf (config, " sx1(1/ns) sy1(1/ns) sz1(1/ns)\n"); for (i=0; i<3*np; i+=3) fprintf (config, "%s %f %f %f %f %.6f %.6f %.6f\n", tp[i/3], mass[i/3]* (39.948*1E-3/ 6.0221367E+23 ) *1000* 6.0221367E+23 , s[i], s[i+1], s[i+2], s1[i] /delta/ (((3.405E-10) *sqrt((39.948*1E-3/ 6.0221367E+23 ) / (119.8* 1.380658E-23 ) )) *1E12) *1000, s1[i+1]/delta/ (((3.405E-10) *sqrt((39.948*1E-3/ 6.0221367E+23 ) / (119.8* 1.380658E-23 ) )) *1E12) *1000, s1[i+2]/delta/ (((3.405E-10) *sqrt((39.948*1E-3/ 6.0221367E+23 ) / (119.8* 1.380658E-23 ) )) *1E12) *1000); fclose (config); return; } void save_tmp_config () { char buf[128 ]; static int tmp_config_counter = 1; if (strcasecmp(tmp_config, "NULL")) { sprintf (buf, "%s%d", tmp_config, tmp_config_counter); write_config (buf); printf ("saved on \"%s\" at step %d, t = %.2f ps.\n\n", buf, kb, kb*delta* (((3.405E-10) *sqrt((39.948*1E-3/ 6.0221367E+23 ) / (119.8* 1.380658E-23 ) )) *1E12) ); tmp_config_counter++; } return; } void save_mesh (char filename[]) { int i; FILE *out; double tot = volume / 100 * (kb-kstart+1); if ( strcasecmp(filename, "NULL") ) { out = fopen (filename, "w"); fprintf (out, "%d %f %d %d %f %f\n", kb, kb*delta* (((3.405E-10) *sqrt((39.948*1E-3/ 6.0221367E+23 ) / (119.8* 1.380658E-23 ) )) *1E12) , kstart, kb-kstart+1, H[2][2]*nc[3]/nc[2], H[2][2]*df[2]); for (i=0; i< 100 ; i++) fprintf( out, "%f %f %f %f %f %f\n", (i+0.5)/ 100 *H[2][2], density[i] / tot, density[i]>0? vx[i] / density[i] : 0, density[i]>0? (119.8* 1.380658E-23 ) / 1.380658E-23 * (temp[i]/density[i]- (39.948*1E-3/ 6.0221367E+23 ) / (39.948*1E-3/ 6.0221367E+23 ) * (( vx[i]/density[i] )*( vx[i]/density[i] )) )/3. : 0, density[i]>0? taum[i][0][2]/tot : 0, density[i]>0? (taum[i][0][0]+taum[i][1][1]+ taum[i][2][2])/3./tot : 0 ); } fclose(out); return; } void ds (double s1[], double s2[], double d[]) { d[0] = s2[0] - s1[0]; d[1] = s2[1] - s1[1]; d[2] = s2[2] - s1[2]; while (d[0] > 0.5) d[0]--; while (d[0] < -0.5) d[0]++; while (d[1] > 0.5) d[1]--; while (d[1] < -0.5) d[1]++; while (d[2] > 0.5) d[2]--; while (d[2] < -0.5) d[2]++; return; } double r2 (int i, int j, double si[], double sj[]) { double ds1, ds2, ds3; ds1 = sj[3*j] - si[3*i]; ds2 = sj[3*j+1] - si[3*i+1]; ds3 = sj[3*j+2] - si[3*i+2]; if (ds1 > 0.5) ds1--; else if (ds1 < -0.5) ds1++; if (ds2 > 0.5) ds2--; else if (ds2 < -0.5) ds2++; if (ds3 > 0.5) ds3--; else if (ds3 < -0.5) ds3++; r[0] = ds1*H[0][0] + ds2*H[1][0] + ds3*H[2][0]; r[1] = ds1*H[0][1] + ds2*H[1][1] + ds3*H[2][1]; r[2] = ds1*H[0][2] + ds2*H[1][2] + ds3*H[2][2]; return (r[3] = r[0]*r[0] + r[1]*r[1] + r[2]*r[2]); } void pair_potential() { int i, j, k, l; double r1, ir2, ir6, ir12, a, b, vij; double ff, ffx, ffy, ffz, cc, dd; double res = LJ6_12?(4*(pow((2.5* (3.405E-10) / (3.405E-10) ) ,-12.)-pow((2.5* (3.405E-10) / (3.405E-10) ) ,-6.))) : ((0.4228* 4.3597482E-18 /1000./ (119.8* 1.380658E-23 ) ) *(0.6060 *exp(- sqrt((0.6060 +1)* (0.6627* 4.3597482E-18 /1000./ (119.8* 1.380658E-23 ) / pow(0.529177249 / ((3.405E-10) *1E10) , 2.)) / 0.6060 / (0.4228* 4.3597482E-18 /1000./ (119.8* 1.380658E-23 ) ) ) *((2.5* (3.405E-10) / (3.405E-10) ) - (7.1714* 0.529177249 / ((3.405E-10) *1E10) ) ))- (0.6060 +1)*exp(- sqrt(0.6060 * (0.6627* 4.3597482E-18 /1000./ (119.8* 1.380658E-23 ) / pow(0.529177249 / ((3.405E-10) *1E10) , 2.)) /(0.6060 +1)/ (0.4228* 4.3597482E-18 /1000./ (119.8* 1.380658E-23 ) ) ) *((2.5* (3.405E-10) / (3.405E-10) ) - (7.1714* 0.529177249 / ((3.405E-10) *1E10) ) )))) ; double ref = LJ6_12?((48*pow((2.5* (3.405E-10) / (3.405E-10) ) ,-12.)-24*pow((2.5* (3.405E-10) / (3.405E-10) ) ,-6.))/ (2.5* (3.405E-10) / (3.405E-10) ) ) : ((0.4228* 4.3597482E-18 /1000./ (119.8* 1.380658E-23 ) ) *(0.6060 * sqrt((0.6060 +1)* (0.6627* 4.3597482E-18 /1000./ (119.8* 1.380658E-23 ) / pow(0.529177249 / ((3.405E-10) *1E10) , 2.)) / 0.6060 / (0.4228* 4.3597482E-18 /1000./ (119.8* 1.380658E-23 ) ) ) *exp(- sqrt((0.6060 +1)* (0.6627* 4.3597482E-18 /1000./ (119.8* 1.380658E-23 ) / pow(0.529177249 / ((3.405E-10) *1E10) , 2.)) / 0.6060 / (0.4228* 4.3597482E-18 /1000./ (119.8* 1.380658E-23 ) ) ) *((2.5* (3.405E-10) / (3.405E-10) ) - (7.1714* 0.529177249 / ((3.405E-10) *1E10) ) ))- (0.6060 +1)* sqrt(0.6060 * (0.6627* 4.3597482E-18 /1000./ (119.8* 1.380658E-23 ) / pow(0.529177249 / ((3.405E-10) *1E10) , 2.)) /(0.6060 +1)/ (0.4228* 4.3597482E-18 /1000./ (119.8* 1.380658E-23 ) ) ) *exp(- sqrt(0.6060 * (0.6627* 4.3597482E-18 /1000./ (119.8* 1.380658E-23 ) / pow(0.529177249 / ((3.405E-10) *1E10) , 2.)) /(0.6060 +1)/ (0.4228* 4.3597482E-18 /1000./ (119.8* 1.380658E-23 ) ) ) *((2.5* (3.405E-10) / (3.405E-10) ) - (7.1714* 0.529177249 / ((3.405E-10) *1E10) ) )))) ; pote = kine = cj[0] = cj[1] = cj[2] = 0; for (i=0; i=kstart) ) { j = floor(s[i+2]* 100 ); while (j > 100 ) j -= 100 ; while (j < 0) j += 100 ; density[j] += 1; vx[j] += v[i]; temp[j] += 2 * cc; for (k=0; k<3; k++) for (l=k; l<3; l++) taum[j][k][l] += mass[i/3] * v[i+k] * v[i+l]; } } for (i=0; i 0.5) dd--; l = floor((s[3*i+2]+dd/2)* 100 ); while (l > 100 ) l -= 100 ; while (l < 0) l += 100 ; if ( mb[i] && (!tpd[i]) && mb[j] && (!tpd[j]) && (kb>=kstart) ) { taum[l][0][0] += ffx * r[0]; taum[l][1][1] += ffy * r[1]; taum[l][2][2] += ffz * r[2]; taum[l][0][1] += ffx * r[1]; taum[l][0][2] += ffx * r[2]; taum[l][1][2] += ffy * r[2]; } dd = ( ffx * (v[3*i] + v[3*j] ) + ffy * (v[3*i+1] + v[3*j+1]) + ffz * (v[3*i+2] + v[3*j+2]) ) / 2.; cj[0] += dd * r[0]; cj[1] += dd * r[1]; cj[2] += dd * r[2]; c[0][0] += cc*r[0]*r[0]*r[0]*r[0]; c[0][1] += cc*r[0]*r[0]*r[1]*r[1]; c[0][2] += cc*r[0]*r[0]*r[2]*r[2]; c[1][1] += cc*r[1]*r[1]*r[1]*r[1]; c[1][2] += cc*r[1]*r[1]*r[2]*r[2]; c[2][2] += cc*r[2]*r[2]*r[2]*r[2]; c[0][3] += cc*r[0]*r[0]*r[1]*r[2]; c[0][4] += cc*r[0]*r[0]*r[0]*r[2]; c[0][5] += cc*r[0]*r[0]*r[0]*r[1]; c[1][3] += cc*r[1]*r[1]*r[1]*r[2]; c[1][4] += cc*r[1]*r[1]*r[0]*r[2]; c[1][5] += cc*r[1]*r[1]*r[0]*r[1]; c[2][3] += cc*r[2]*r[2]*r[1]*r[2]; c[2][4] += cc*r[2]*r[2]*r[0]*r[2]; c[2][5] += cc*r[2]*r[2]*r[0]*r[1]; c[3][3] += cc*r[1]*r[2]*r[1]*r[2]; c[3][4] += cc*r[1]*r[2]*r[0]*r[2]; c[3][5] += cc*r[1]*r[2]*r[0]*r[1]; c[4][4] += cc*r[0]*r[2]*r[0]*r[2]; c[4][5] += cc*r[0]*r[2]*r[0]*r[1]; c[5][5] += cc*r[0]*r[1]*r[0]*r[1]; } for (i=0; i 1) s[3*i]--; while (s[3*i+1] < 0) s[3*i+1]++; while (s[3*i+1] > 1) s[3*i+1]--; while (s[3*i+2] < 0) s[3*i+2]++; while (s[3*i+2] > 1) s[3*i+2]--; l = mybin[i] = ((( (int)floor(s[3*i] * maxbin[0]) +(maxbin[0]<<8))%maxbin[0]*maxbin[1]+ ( (int)floor(s[3*i+1] * maxbin[1]) +(maxbin[1]<<8))%maxbin[1])*maxbin[2]+ ( (int)floor(s[3*i+2] * maxbin[2]) +(maxbin[2]<<8))%maxbin[2]) ; safe_append (bindx, binlist, i, l, 0, maxbin[3]); sa[3*i] = s[3*i]; sa[3*i+1] = s[3*i+1]; sa[3*i+2] = s[3*i+2]; } for (i=0; i m ))||((!(( i + m )%2))&&( i < m ))) ) if (tag_top >= 256 ) { printf ("re_anchor: MAX_TAG = %d exceeded.\n", 256 ); exit(1); } else tag[tag_top++] = m; else { for (n=idx[2*m]; (n reference state\n", H[1][0], H[1][1], H[1][2]); printf (" | %f %f %f |\n", H[2][0], H[2][1], H[2][2]); printf ("is divided into %d x %d x %d = %d bins, of thicknesses\n", maxbin[0], maxbin[1], maxbin[2], maxbin[3]); printf ("(%.3f %.3f %.3f) * sqrt(1-2*%.3f) >= RLIST = %.3f.\n\n", dd[0], dd[1], dd[2], MAX_STRAIN, RLIST); particles_move_in(); return; } void print_list_statistics (FILE *out) { int i, j, l, m, n; for (i=j=m=0,n=2L<<20; i bindx[2*i+1]-bindx[2*i]) n = bindx[2*i+1]-bindx[2*i]; j += pow(bindx[2*i+1]-bindx[2*i], 2.); } fprintf (out, "Bin-particle list (%.3f Mbytes allocated, %.1f per " "bin):\n", (double)bindx[2*maxbin[3]]*sizeof(int)/(2L<<20), (double)bindx[2*maxbin[3]]/maxbin[3]); fprintf (out, "max = %d, min = %d, avg occup. = %.1f = %.1f%%, " "fluc = %.2f.\n", m, n, (double)np/maxbin[3], 100.*np/bindx[2*maxbin[3]], sqrt((double)j/maxbin[3]-pow((double)np/maxbin[3],2.))); for (i=j=l=m=0,n=2L<<20; i idx[2*i+1]-idx[2*i]) n = idx[2*i+1]-idx[2*i]; } fprintf (out, "Particle-particle list (%.3f Mbytes allocated, " "%.1f per particle):\n", (double)idx[2*np]*sizeof(int)/(2L<<20), (double)idx[2*np]/np); fprintf (out, "max = %d, min = %d, avg occup. = %.1f = %.1f%%, " "fluc = %.2f.\n\n", m, n, (double)l/np, 100.*l/idx[2*np], sqrt((double)j/np-pow((double)l/np,2.))); return; } void accumulate_and_save_gr (char fn_gr[]) { static int nw_gr = 0; static long igr[200 +1][3] = {0}; static double TRsum = 0, densitysum = 0; int i, j, k, l, m, n, w, gb[3], over, overbin[125 ]; double cc, dd, ee, ff, gg; FILE *gr; for (i=0; i<3; i++) gb[i] = ceil(GRCUT*sqrt(HI[0][i]*HI[0][i]+ HI[1][i]*HI[1][i]+ HI[2][i]*HI[2][i])*maxbin[i]); for (l=0; l= 125 ) { printf ("accumulate_and_save_gr: " "GRMAXBINS = %d exceeded.\n", 125 ); exit(1); } for (w=bindx[2*l]; w binlist[n] ))||((!(( binlist[w] + binlist[n] )%2))&&( binlist[w] < binlist[n] ))) ) if ( r2(binlist[w], binlist[n], s, s) < GRCUT*GRCUT ) ++igr[(int)floor(r[3]/r2del)] [tpd[binlist[w]]+tpd[binlist[n]]]; overbin[over++] = m; } } } nw_gr++; TRsum += T_md_to_T_real(kine* (119.8* 1.380658E-23 ) /1.5/neff/ 1.380658E-23 ); densitysum += np/ 6.0221367E+23 /volume/ pow((3.405E-10) ,3.) ; gr = fopen (fn_gr, "w+"); fprintf (gr, "clf; nw_gr = %d; GR = [\n", nw_gr); for (i=0; i< 200 ; i++) { cc = 4* 3.14159265358979323846 /3*(pow((i+1)*r2del,1.5)-pow(i*r2del,1.5)); dd = cc*nw_gr*(np-nar)*(np-nar-1)/volume/2.; ee = cc*nw_gr*nar*(np-nar)/volume; ff = cc*nw_gr*nar*(nar-1)/volume/2.; gg = dd + ee + ff; if (fabs(dd) > (1.E-8) ) dd = 1./dd; if (fabs(ee) > (1.E-8) ) ee = 1./ee; if (fabs(ff) > (1.E-8) ) ff = 1./ff; if (fabs(gg) > (1.E-8) ) gg = 1./gg; fprintf (gr, "%f %f", cbrt((pow((i+1)*r2del,1.5)+pow(i*r2del,1.5))/2.)* ((3.405E-10) *1E10) , igr[i][0]*dd); if (ntpd > 0) fprintf (gr, " %f %f %f", igr[i][1]*ee, igr[i][2]*ff, (igr[i][0]+igr[i][1]+igr[i][2])*gg); fprintf (gr, "\n"); } fprintf (gr, "]; plot(GR(:,1),GR(:,2));\n"); fprintf (gr, "v = axis; axis([0 %f v(3) v(4)]);\n", GRCUT* ((3.405E-10) *1E10) ); fprintf (gr, "xlabel('r [Angstrom]'); ylabel('g(r)');\n"); fprintf (gr, "title('RDF of Argon at T = %f K, \\rho = %f " "mol/m^3');\n", TRsum/nw_gr, densitysum/nw_gr); if (nar > 0) { fprintf (gr, "hold on; plot(GR(:,1),GR(:,3),'r');\n"); fprintf (gr, "plot(GR(:,1),GR(:,4),'g');\n"); fprintf (gr, "plot(GR(:,1),GR(:,5),'b');\n"); } fclose (gr); return; } void print_atom_statistics (FILE *out, char *sp) { int j; fprintf(out,"%sThere are %d atoms other than \"%2s of %.4f amu\".\n", sp, nk, "Ar" , (39.948*1E-3/ 6.0221367E+23 ) *1000* 6.0221367E+23 ); fprintf(out,"%saltogether %d different kinds will coexist:\n", sp, nkind); fprintf(out,"%s--------------------------------------\n", sp); fprintf(out,"%s Type Mass(amu) Occurrences\n", sp); for (j=0; j 0) { fprintf(out,"%samong which %d kinds: ", sp, ntpd); for (j=0; j( l_ratio[1] )?( l_ratio[0] ):( l_ratio[1] )) ; maxl_ratio = (( maxl_ratio )>( l_ratio[2] )?( maxl_ratio ):( l_ratio[2] )) ; fprintf(out,"%sBecause of it, the system bin size has to increase\n" "%sfrom pure-%2s by maxl_ratio = %f; and\n", sp, sp, "Ar" , maxl_ratio); minl_ratio[0] = (( l_ratio[0] )<( l_ratio[1] )?( l_ratio[0] ):( l_ratio[1] )) ; minl_ratio[1] = (( l_ratio[2] )<( l_ratio[1] )?( l_ratio[2] ):( l_ratio[1] )) ; fprintf(out,"%sMAX_DRIFT's of %2s, N%2s respectively should decrease " "by\n", sp, "Ar" , "Ar" ); fprintf(out,"%sminl_ratio[0] = %f, minl_ratio[1] = %f.\n", sp, minl_ratio[0], minl_ratio[1]); } if (nmb != np) { fprintf(out,"\n%sAnd, we have found %d atoms with nominal " "masses\n", sp, np-nmb); fprintf(out,"%s> %.3f AMU, and they will be treated as " "immobile.\n", sp, (9999.) ); } fprintf(out,"\n"); waitfor ((0.) ); return; } int TFE (double sz1, double sz2) { int i, j, k; double cc, vxavg, kei, sz[2* ((( 1 )>( 2 )?( 1 ):( 2 )) ) +1], vx_sz[1 +1], ke_sz[2 +1], S[2 +1][2 +1], SI[2 +1][2 +1]; for (i=0; i<2* 1 +1; i++) sz[i] = 0.; for (i=0; i< 1 +1; i++) vx_sz[i] = 0.; for (i=0; i= (( sz1 )<( sz2 )?( sz1 ):( sz2 )) ) (s[i+2] <= (( sz1 )>( sz2 )?( sz1 ):( sz2 )) ) ) for (cc=1,k=0; k<2* 1 +1; k++,cc*=s[i+2]) { sz[k] += cc; if (k< 1 +1) vx_sz[k] += v[i] * cc; } if (sz[0] > 0) { sz[1] /= sz[0]; sz[2] /= sz[0]; vx_sz[0] /= sz[0]; vx_sz[1] /= sz[0]; af[0] = ( vx_sz[0]*sz[2] - sz[1]*vx_sz[1] ) / (sz[2] - sz[1]*sz[1]); af[1] = ( vx_sz[1] - sz[1]*vx_sz[0] ) / (sz[2] - sz[1]*sz[1]); for (i=0; i<2* 2 +1; i++) sz[i] = 0.; for (i=0; i< 2 +1; i++) ke_sz[i] = 0.; for (i=0; i= (( sz1 )<( sz2 )?( sz1 ):( sz2 )) ) (s[i+2] <= (( sz1 )>( sz2 )?( sz1 ):( sz2 )) ) ) { vxavg = af[0] + af[1] * s[i+2]; kei = mass[i/3] / 3 * (119.8* 1.380658E-23 ) / 1.380658E-23 * ( (( v[i+1] )*( v[i+1] )) + (( v[i+2] )*( v[i+2] )) + (( v[i]-vxavg )*( v[i]-vxavg )) ); for (cc=1,k=0; k<2* 2 +1; k++,cc*=s[i+2]) { sz[k] += cc; if (k< 2 +1) ke_sz[k] += kei * cc; } for (j=0; j< 2 +1; j++) for (k=0; k< 2 +1; k++) S[j][k] = s[j+k]; matinv (S, SI); for (j=0; j< 2 +1; j++) for (af[2+j]=k=0; k< 2 +1; k++) af[2+j] += SI[j][k] * ke_sz[k]; } } return ((int)sz[0]); } int main (int argc, char *argv[]) { int i, j, k, l, m, n; double cc, dd, ee, ff, tmp, *calm; struct schedule_list *sch; for (buf[0]=(char)0,i=0; i= maxkb) { printf ("error: kstart >= MAXKB, " "leaving no room for averaging.\n"); exit(1); } else if (kstart >= 2*maxkb/3) { printf ("** Warning: avg. period small compared with MAXKB.\n\n"); waitfor ((0.) ); } fscanf (((_IO_FILE*)(&_IO_stdin_)) , "Load configuration from file ('NULL'->start anew):\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%s\n\n", fn_config_read); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "NTH or NTS mode, VOLUMETRIC or free:\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%s", buf); NTH = !strcasecmp(buf, "NTH"); NTS = !strcasecmp(buf, "NTS"); if ( (!NTH) && (!NTS) ) { printf("Error: \"%s\" run mode is invalid,\n", buf); printf("must select between NTH or NTS mode.\n"); return(1); } fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%s\n\n", buf); VOLUMETRIC = !strcasecmp(buf, "VOLUMETRIC"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "Default schedule of the reference unit cell [A]:\n"); for (HSCH=HCALM=i=0; i<3; i++) for (j=i; j<3; j++) { for (k=0; (buf[k]= _IO_getc ( ((_IO_FILE*)(&_IO_stdin_)) ) )!='\n'; k++); buf[k] = '\0'; sch = add_schedule (Hschname[((5- (( i )<( j )?( i ):( j )) )* (( i )<( j )?( i ):( j )) /2+ (( i )>( j )?( i ):( j )) ) ], buf, 0, maxkb); HSCH = HSCH || sch->nontrivial; calm = calm_period (Hschname[((5- (( i )<( j )?( i ):( j )) )* (( i )<( j )?( i ):( j )) /2+ (( i )>( j )?( i ):( j )) ) ], kstart, maxkb); HCALM = (calm== ((void *)0) )? maxkb : (( HCALM )>( calm[0] )?( HCALM ):( calm[0] )) ; } if (NTS) HSCH = 0; fscanf (((_IO_FILE*)(&_IO_stdin_)) , "\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "nc(1,2,3), number of wall layers among nc(3):\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%d %d %d %d\n\n", &nc[0], &nc[1], &nc[2], &nc[3]); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "OPC vx [reduced], T [K], sz_OPC, sz_TFE_bottom:\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%lf %lf %lf %lf\n\n", df[0], df[1], df[2], df[3]); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "Default multiplications of the unit cell " "(also structure factor):\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%d %d %d\n\n", &nc[0], &nc[1], &nc[2]); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "In NTH, use the above H schedule " "instead of the configuration file:\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%s\n\n", buf); if ( strcasecmp(fn_config_read, "NULL") ) { USE_CONFIG_H = NTS || ((buf[0]!='1') && (buf[0]!='y') && (buf[0]!='Y')); if ( USE_CONFIG_H ) { HSCH = 0; HCALM = 0; printf ("Do NOT override H from \"%s\" --\n" "default H schedule abandoned.\n\n", fn_config_read); } else printf ("OVERRIDE H from \"%s\" --\n" "will use default H schedule (%s).\n\n", fn_config_read, HSCH?(HCALM<=kstart?"nontrivial but calm": "nontrivial and not calm"):"trivial"); waitfor ((0.) ); } fscanf (((_IO_FILE*)(&_IO_stdin_)) , "Default external stress [MPa] schedule:\n"); for (SSCH=SCALM=i=0; i<3; i++) for (j=i; j<3; j++) { for (k=0; (buf[k]= _IO_getc ( ((_IO_FILE*)(&_IO_stdin_)) ) )!='\n'; k++); buf[k] = '\0'; sch = add_schedule (Sschname[((5- (( i )<( j )?( i ):( j )) )* (( i )<( j )?( i ):( j )) /2+ (( i )>( j )?( i ):( j )) ) ], buf, 0, maxkb); SSCH = SSCH || sch->nontrivial; calm = calm_period (Sschname[((5- (( i )<( j )?( i ):( j )) )* (( i )<( j )?( i ):( j )) /2+ (( i )>( j )?( i ):( j )) ) ], kstart, maxkb); SCALM = (calm== ((void *)0) )? maxkb : (( SCALM )>( calm[0] )?( SCALM ):( calm[0] )) ; } if (NTH) SSCH = 0; fscanf (((_IO_FILE*)(&_IO_stdin_)) , "\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "Save configuration on file:\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%s\n\n", fn_config_write); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "Pair potential to be used ('LJ6_12' or 'WOON'):\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%s\n\n", potfun_name); LJ6_12 = !strcasecmp (potfun_name, "LJ6_12"); WOON = !strcasecmp (potfun_name, "WOON"); if ( (!LJ6_12) && (!WOON) ) { printf("Error: select between LJ6_12 or WOON potential.\n"); exit(1); } fscanf (((_IO_FILE*)(&_IO_stdin_)) , "Request calculating thermal conductivity " "and shear viscosity?\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%s\n\n", buf); CALC = (buf[0]=='1') || (buf[0]=='y') || (buf[0]=='Y'); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "MAX_STRAIN for strain-triggered " "bin re-partition:\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%s\n\n", buf); MAX_STRAIN = strcasecmp(buf, "default")? fabs(atof(buf)) : (NTH&&(!HSCH))? 0.001 : 0.02 ; fscanf (((_IO_FILE*)(&_IO_stdin_)) , "MAX_DRIFT for particle re-anchoring:\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%s\n\n", buf); MAX_DRIFT = strcasecmp(buf, "default")?fabs(atof(buf)): 0.08 ; RLIST = sqrt((1+2*MAX_STRAIN)/(1-2*MAX_STRAIN)) * (1+2*MAX_DRIFT)* (2.5* (3.405E-10) / (3.405E-10) ) ; fscanf (((_IO_FILE*)(&_IO_stdin_)) , "GRCUT for g(r) calculation [A]:\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%s\n\n", buf); GRCUT = strcasecmp(buf, "default")? fabs(atof(buf)): (15.) ; GRCUT /= ((3.405E-10) *1E10) ; fscanf (((_IO_FILE*)(&_IO_stdin_)) , "Temperature [K] schedule:\n"); for (k=0; (buf[k]= _IO_getc ( ((_IO_FILE*)(&_IO_stdin_)) ) )!='\n'; k++); buf[k] = '\0'; sch = add_schedule ("REALT", buf, 0, maxkb); TSCH = sch->nontrivial; calm = calm_period ("REALT", kstart, maxkb); TCALM = (calm== ((void *)0) )? maxkb : calm[0]; fscanf (((_IO_FILE*)(&_IO_stdin_)) , "\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "Integration timestep size [fs]:\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%lf\n\n", &delta); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "Temperature/stress coupling time constant [ps]:\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%s\n\n", buf); tau = ( strncasecmp(buf,"inf",3) && strcasecmp(buf,"huge") )? fabs(atof(buf) ): -1; fscanf (((_IO_FILE*)(&_IO_stdin_)) , "Frequency of output on screen and disk:\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%d\n\n", &kw_lp); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "Random number generator seed (0->time seed):\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%d\n\n", &iseed); if (iseed == 0) iseed = time(((void *)0) ); srandom ((unsigned)iseed); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "Create voids at (sx sy sz radius[A]):\n"); for (nvoids=0; ; nvoids++) { while ( ((buf[0]= _IO_getc ( ((_IO_FILE*)(&_IO_stdin_)) ) ) != '\n') && (buf[0] != (char)(-1) ) && (buf[0] != 'n') && (buf[0] != 'N') && (buf[0] != '(') ); if (buf[0] == '(') { if (nvoids >= 128 ) { printf ("Error: MAX_VOIDS = %d exceeded.\n", 128 ); exit(1); } fscanf (((_IO_FILE*)(&_IO_stdin_)) , "%lf %lf %lf %lf)", &voids[nvoids][0], &voids[nvoids][1], &voids[nvoids][2], &voids[nvoids][3]); voids[nvoids][3] /= ((3.405E-10) *1E10) ; } else { if ((buf[0] == 'n')||(buf[0] == 'N')) while ( ((buf[0]= _IO_getc ( ((_IO_FILE*)(&_IO_stdin_)) ) != '\n')) && (buf[0] != (char)(-1) ) ); break; } } fscanf (((_IO_FILE*)(&_IO_stdin_)) , "\n"); fscanf (((_IO_FILE*)(&_IO_stdin_)) , "Save on temporary config (NULL=don't) " "with frequency:\n%s %s\n", tmp_config, buf); if (!strcasecmp(tmp_config, "default")) strcpy(tmp_config, "tmp." ); tmp_config_freq = (!strcasecmp(buf, "default"))? 10 *kw_lp: atoi(buf); printf ("Instructions read complete.\n\n"); if (strcasecmp(fn_config_read, "NULL")) { if ((config=fopen(fn_config_read,"r")) == ((void *)0) ) { printf("Error: cannot open config file \"%s\".\n", fn_config_read); exit(1); } printf ("The initial configuration is loaded from \"%s\".\n", fn_config_read); fscanf (config, "Number of particles = %d\n\n", &np); } else { printf ("Initial config. is structurally perfect " "fcc crystal.\n"); np = (6) *nc[0]*nc[1]*nc[2]; } printf ("Number of particles in this simulation = %d.\n\n", np); np3 = 3 * np; sa = (double *) malloc (np3*sizeof(double)); s = (double *) malloc (np3*sizeof(double)); s1 = (double *) malloc (np3*sizeof(double)); s2 = (double *) malloc (np3*sizeof(double)); s3 = (double *) malloc (np3*sizeof(double)); s4 = (double *) malloc (np3*sizeof(double)); s5 = (double *) malloc (np3*sizeof(double)); v = (double *) malloc (np3*sizeof(double)); f = (double *) malloc (np3*sizeof(double)); fs = (double *) malloc (np3*sizeof(double)); d = (double *) malloc (np3*sizeof(double)); mass = (double *) malloc (np*sizeof(double)); ep = (double *) malloc (np*sizeof(double)); mybin = (int *) malloc (np*sizeof(int)); tp = (char **) malloc (np*sizeof(char *)); tp[0] = (char *) malloc (np*(2 +1)*sizeof(char)); for (i=1; i=0; j--) { if ( USE_CONFIG_H && (j>=i) ) { cancel_schedule (Hschname[((5- (( i )<( j )?( i ):( j )) )* (( i )<( j )?( i ):( j )) /2+ (( i )>( j )?( i ):( j )) ) ]); sprintf (buf, "%f", (H[i][j]/nc[i]+H[j][i]/nc[j])/2); add_schedule (Hschname[((5- (( i )<( j )?( i ):( j )) )* (( i )<( j )?( i ):( j )) /2+ (( i )>( j )?( i ):( j )) ) ], buf, 0, maxkb); } H[i][j] = nc[i] * schedule(Hschname[((5- (( i )<( j )?( i ):( j )) )* (( i )<( j )?( i ):( j )) /2+ (( i )>( j )?( i ):( j )) ) ], 0) / ((3.405E-10) *1E10) ; H1[i][j] = H2[i][j] = H3[i][j] = H4[i][j] = H5[i][j] = 0.; } volume = matinv (H, HI); mateqv (HI, HIBEGIN); if (HSCH) printf("There is a nontrivial H schedule, with\n"); printf (" | %f %f %f |\n", H[0][0]* ((3.405E-10) *1E10) , H[0][1]* ((3.405E-10) *1E10) , H[0][2]* ((3.405E-10) *1E10) ); printf ("Initial H = | %f %f %f | A, nc = (%d %d %d)\n", H[1][0]* ((3.405E-10) *1E10) , H[1][1]* ((3.405E-10) *1E10) , H[1][2]* ((3.405E-10) *1E10) , nc[0], nc[1], nc[2]); printf (" | %f %f %f |\n", H[2][0]* ((3.405E-10) *1E10) , H[2][1]* ((3.405E-10) *1E10) , H[2][2]* ((3.405E-10) *1E10) ); if (nc[0]*nc[1]*nc[2] <= 0) { printf ("Error: illegal nc[] partition of H into " "reference \"unit cells\"."); exit(1); } if ( HSCH && (HCALM <= kstart) ) printf("but which stops changing after KSTART;\n"); for (i=k=0; i k ) { printf ("\nTo create %d voids, we throw away %d particles,\n", nvoids, np-k); printf ("and now has only %d particles.\n\n", np=k); waitfor ((0.) ); } l_ratio[0] = 1.; l_ratio[1] = sqrt((1.0) ) ; l_ratio[2] = (1.0) ; e_ratio[0] = 1.; e_ratio[1] = (((1.0) +1)/2.) ; e_ratio[2] = (1.0) ; for (nkind=i=0; i= (9999.) - (1.E-8) ) s1[3*i] = s1[3*i+1] = s1[3*i+2] = mb[i] = 0; else {nmb++; mb[i] = 1;} neff = nmb - 1; maxl_ratio = minl_ratio[0] = minl_ratio[1] = 1.; if (nkind > 1) print_atom_statistics(((_IO_FILE*)(&_IO_stdout_)) , ""); printf ("Initial reduced density = %f = %.1f mol/m^3.\n\n", np/volume, np/volume/ 6.0221367E+23 / pow((3.405E-10) ,3.) ); waitfor ((0.) ); if (NTH) { printf ("Simulation will be carried out in NTH mode\n"); if (VOLUMETRIC) printf ("keeping the initial volume invariant;\n"); } else { printf ("Simulation will be carried out in NTS mode\n"); for (i=0; i<3; i++) for (j=i; j<3; j++) { cc = schedule(Sschname[((5- (( i )<( j )?( i ):( j )) )* (( i )<( j )?( i ):( j )) /2+ (( i )>( j )?( i ):( j )) ) ],0) / (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ; sout[i][j] = sout[j][i] = cc; } if (VOLUMETRIC) { printf ("with hydrostatic constraints;\n"); trace_decompose (sout, sout, ((void *)0) ); } if (SSCH) printf("There is a nontrivial stress schedule with initial\n"); else printf("The loading stay constant during the simulation, with\n"); printf (" | %f %f %f |\n", sout[0][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , sout[0][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , sout[0][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); printf ("external stress = | %f %f %f | MPa,\n", sout[1][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , sout[1][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , sout[1][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); printf (" | %f %f %f |\n", sout[2][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , sout[2][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , sout[2][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); if ( SSCH && (SCALM <= kstart) ) printf("but which stops changing after KSTART;\n"); hydro = trace_decompose (sout, ((void *)0) , sout); shear = sqrt( sout[0][0]*sout[0][0]/2 + sout[1][1]*sout[1][1]/2 + sout[2][2]*sout[2][2]/2 + sout[0][1]*sout[0][1] + sout[0][2]*sout[0][2] + sout[1][2]*sout[1][2] ); printf ("mean normal stress = %f MPa,\n", hydro* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); printf ("von Mises shear stress = %f MPa,\n", shear* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); } printf ("using the %s pair-potential.\n\n", potfun_name); waitfor ((0.) ); if (TSCH) printf("We are planning for a nontrivial T schedule with\n"); TR = schedule ("REALT", 0); printf ("Initial real temperature = %.2f K\n", TR); T = T_real_to_T_md (TR, &dTdR); if (TSCH) printf ("and final real temperature = %.2f K;\n", schedule ("REALT", maxkb)); if ( TSCH && (TCALM <= kstart) ) printf("but the schedule stops changing after KSTART.\n"); printf("\n"); waitfor ((0.) ); wallmass = (tau>0)? (2000.) / (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) * cbrt(volume) * (( tau/ (((3.405E-10) *sqrt((39.948*1E-3/ 6.0221367E+23 ) / (119.8* 1.380658E-23 ) )) *1E12) )*( tau/ (((3.405E-10) *sqrt((39.948*1E-3/ 6.0221367E+23 ) / (119.8* 1.380658E-23 ) )) *1E12) )) : 1000. * (39.948*1E-3/ 6.0221367E+23 ) / (39.948*1E-3/ 6.0221367E+23 ) ; printf ("MAXKB = %d\t KSTART = %d\n", maxkb, kstart); printf ("KW_LP = %d\t ISEED = %d\n", kw_lp, iseed); printf ("delta = %.2f ns tau = %.2f ps\n", delta, tau); printf ("wallmass = %.2f Argon mass\n\n", wallmass * (39.948*1E-3/ 6.0221367E+23 ) / (39.948*1E-3/ 6.0221367E+23 ) ); waitfor ((0.) ); delta /= ((((3.405E-10) *sqrt((39.948*1E-3/ 6.0221367E+23 ) / (119.8* 1.380658E-23 ) )) *1E12) *1000); tau /= (((3.405E-10) *sqrt((39.948*1E-3/ 6.0221367E+23 ) / (119.8* 1.380658E-23 ) )) *1E12) ; printf ("Across the entire simulation,\n"); if (tau < 0) if (NTH) if (!HSCH) printf ("internal energy is conserved because the\n" "temperature coupling is turned off and\n" "there is no active H schedule.\n\n"); else printf ("internal energy may not be conserved although\n" "the temperature coupling is turned off because\n" "there is an active H schedule.\n\n"); else if (!SSCH) if (VOLUMETRIC) printf ("the hydrostatic enthalpy is conserved because\n" "the temperature coupling is turned off, there\n" "is NO active stress schedule + VOLUMETRIC.\n\n"); else printf ("generalized enthalpy is conserved (piecewise)\n" "because the temperature coupling is turned off\n" "and there is no active stress schedule.\n\n"); else if (VOLUMETRIC) printf ("the hydrostatic enthalpy may not be conserved\n" "although the temperature coupling is turned off\n" "because there is an active stress schedule.\n\n"); else printf ("the generalized enthalpy may not be conserved\n" "although the temperature coupling is turned off\n" "because there is an active stress schedule.\n\n"); else printf ("no dynamical quantities are conserved because\n" "the temperature coupling is turned on.\n\n"); waitfor ((0.) ); printf ("After KSTART = %d,\n", kstart); if (NTH) if ( (TCALM <= kstart) && (HCALM <= kstart) ) printf ("the conserved quantity is the internal energy\n" "because there are no further T or H activities\n" "and the system automatically switch to NEH mode.\n\n"); else printf ("there are no conserved dynamical quantities\n" "because there are still activities in %s.\n\n", (TCALM<=kstart)?"H":(HCALM<=kstart)?"T":"T and H"); else if ( (TCALM <= kstart) && (SCALM <= kstart) ) printf ("the conserved quantity is the internal energy\n" "as the system switches to NEH mode, since there\n" "are no further T or S activities.\n\n"); else printf ("there are no conserved dynamical quantities\n" "because there are still activities in %s.\n\n", (TCALM<=kstart)?"S":(SCALM<=kstart)?"T":"T and S"); waitfor ((0.) ); if ( (!(NTH && (TCALM <= kstart) && (HCALM <= kstart))) && (!(NTS && (TCALM <= kstart) && (SCALM <= kstart))) ) { printf("************************ WARNING! ************************\n"); printf("Because the system cannot switch to NEH mode after KSTART,\n"); printf("micro-canonical fluctuation formulas does not give the\n"); printf("correct results for second-order thermodynamic derivatives\n"); printf("that appear in the final printouts;\n"); if ( CALC ) { CALC = 0; printf("\nFurthermore, the request for thermal conductivity\n" "calculation is DENIED because it requires NEH modes.\n"); } printf("************************ WARNING! ************************\n"); printf ("\n"); } else if ( CALC ) printf("Thermal conductivity and shear viscosity will be\n" "correctly calculated.\n\n"); waitfor ((0.) ); printf ("Property averages => \"%s\";\n", fn_property); property = fopen (fn_property, "w+"); printf ("Structure evolution => \"%s\";\n", fn_structure); structure = fopen (fn_structure, "w+"); printf ("Temporary configur. => \"%s*\" (every %d steps);\n", tmp_config, tmp_config_freq); printf ("Radial distribution => \"%s\";\n", "gr.m" ); if ( CALC ) { printf ("Heat current => \"%s\", \"%s\", \"%s\";\n", fn_flux[0 ][0], fn_flux[0 ][1], fn_flux[0 ][2]); printf ("Shear stress => \"%s\", \"%s\", \"%s\";\n", fn_flux[1 ][0], fn_flux[1 ][1], fn_flux[1 ][2]); cc = delta * (((3.405E-10) *sqrt((39.948*1E-3/ 6.0221367E+23 ) / (119.8* 1.380658E-23 ) )) *1E12) ; for (i=0; i<3; i++) { heat_current[i] = fopen (fn_flux[0 ][i],"w+"); fwrite (&cc, sizeof(double), 1, heat_current[i]); shear_stress[i] = fopen (fn_flux[1 ][i],"w+"); fwrite (&cc, sizeof(double), 1, shear_stress[i]); } } printf ("Final configuration => \"%s\".\n\n", fn_config_write); waitfor ((0.) ); for (cc=totalmass=0.,i=0; i\n", fn_config_read, cc, (0.1) ); printf ("Assign random velocities to particles ...\n"); x[0] = x[1] = x[2] = 0; randnorm (np3, v); for (i=0; i= 0. ) && (df[2] <= 1. ) && (i=TFE(df[3], df[2])) ) { printf ("By doing linear regression on the %d %s atoms " "inside\n", i, "Ar" ); printf ("sz = (%.5f, %.5f) => vx(%.5f) = %f,\n", df[3], df[2], df[2], cc=af[0]+af[1]*df[2]); printf ("vx(%.5f=%.5f A) = 0 (wall at %.5f=%.5f A).\n\n", -af[0]/af[1], -af[0]/af[1]*H[2][2]* (3.405E-10) *1E10, (double)nc[3]/nc[2], (double)nc[3]/nc[2]* (3.405E-10) *1E10); for (i=0; i possible shrinkage ratio because of that = %f;\n", sqrt(1-2*MAX_STRAIN)/sqrt(1+2*MAX_STRAIN)); printf ("=> RLIST = %.3f = %.3f A = %.3f RCUT.\n\n", RLIST, RLIST* ((3.405E-10) *1E10) , RLIST/ (2.5* (3.405E-10) / (3.405E-10) ) ); waitfor ((0.) ); initialize_all_lists(); nw_repartition = nw_anchor = 0; print_list_statistics (((_IO_FILE*)(&_IO_stdout_)) ); waitfor ((0.) ); i = (( maxbin[0] )<( maxbin[1] )?( maxbin[0] ):( maxbin[1] )) ; i = (( i )<( maxbin[2] )?( i ):( maxbin[2] )) ; if (2*GRCUT > i*maxl_ratio*RLIST) GRCUT = i*maxl_ratio*RLIST/2.; cc = np/2.*(np-1)/volume*4* 3.14159265358979323846 *GRCUT*GRCUT*GRCUT/3/ 200 ; dd = ceil(1./ 0.001 / 0.001 /cc); kw_gr = ceil((maxkb-kstart)/dd); r2del = GRCUT*GRCUT/ 200 ; printf ("GRCUT = %.3f A\tGRMESH = %d\n", GRCUT* ((3.405E-10) *1E10) , 200 ); printf ("GRACCU = %.2f%%\t\tkw_gr = %d\n", 100* 0.001 , kw_gr); printf ("nc = [%d %d %d] (structure factor)\n\n", nc[0], nc[1], nc[2]); waitfor ((0.) ); for (i=0; i( TCALM )?( SCALM ):( TCALM )) + 0.5*kstart; for (i=0; i<6; i++) for (j=0; j<6; j++) { ck1tum[i][j] = 0.; ck1sum[i][j] = 0.; } printf ("Simulation starts ...\n"); for (kb=0; kb=kstart))) ) || ( HSCH && (!((HCALM<=kstart)&&(kb>=kstart))) ); scale_T = (kb=kstart))) ); T = T_real_to_T_md (TR=schedule("REALT",kb), &dTdR); if ( NTH && scale_H ) { for (i=0; i<3; i++) for (j=0; j<3; j++) H[i][j] = schedule (Hschname[((5- (( i )<( j )?( i ):( j )) )* (( i )<( j )?( i ):( j )) /2+ (( i )>( j )?( i ):( j )) ) ],kb) * nc[i] / ((3.405E-10) *1E10) ; if (VOLUMETRIC) { cc = matinv (H, HI); for (i=0; i<3; i++) for (j=0; j<3; j++) { H[i][j] = H[i][j] * cbrt(volume/cc); HI[i][j] = HI[i][j] / cbrt(volume/cc); } } else volume = matinv (H, HI); } else if ( NTS && scale_H ) { for (i=0; i<3; i++) for (j=0; j<3; j++) sout[i][j] = schedule (Sschname[((5- (( i )<( j )?( i ):( j )) )* (( i )<( j )?( i ):( j )) /2+ (( i )>( j )?( i ):( j )) ) ],kb) / (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ; if (VOLUMETRIC) hydro = trace_decompose (sout, sout, ((void *)0) ); else hydro = trace_decompose (sout, ((void *)0) , sout); } for (msd=i=0; i (( minl_ratio[tpd[i]]*MAX_DRIFT* (2.5* (3.405E-10) / (3.405E-10) ) )*( minl_ratio[tpd[i]]*MAX_DRIFT* (2.5* (3.405E-10) / (3.405E-10) ) )) ) { d[3*i] += r[0]; d[3*i+1] += r[1]; d[3*i+2] += r[2]; re_anchor(i); nw_anchor++; anchortum++; msd += d[3*i]*d[3*i]+d[3*i+1]*d[3*i+1]+d[3*i+2]*d[3*i+2]; } else msd += (d[3*i]+r[0])*(d[3*i]+r[0])+ (d[3*i+1]+r[1])*(d[3*i+1]+r[1])+ (d[3*i+2]+r[2])*(d[3*i+2]+r[2]); else if (idx[2*i+1]==idx[2*i+2]) re_anchor(i); } msd /= neff; if (kb==kstart) msd_kstart = msd; if ( NTS ) { matmul (HIBEGIN, H, A); matran (A, B); matmul (A, B, G); for (cc=i=0; i<3; i++) for (j=0; j<3; j++) cc += (( G[i][j]-((i==j)?1.:0.) )*( G[i][j]-((i==j)?1.:0.) )) / 4.; if (cc >= (3* (( 0.5 )*( 0.5 )) ) ) { printf ("\nAt step %d, t = %.2f ps\n", kb, kb * delta * (((3.405E-10) *sqrt((39.948*1E-3/ 6.0221367E+23 ) / (119.8* 1.380658E-23 ) )) *1E12) ); printf (" | %f %f %f |\n", H[0][0]* ((3.405E-10) *1E10) , H[0][1]* ((3.405E-10) *1E10) , H[0][2]* ((3.405E-10) *1E10) ); printf ("H = | %f %f %f | A,\n", H[1][0]* ((3.405E-10) *1E10) , H[1][1]* ((3.405E-10) *1E10) , H[1][2]* ((3.405E-10) *1E10) ); printf (" | %f %f %f |\n", H[2][0]* ((3.405E-10) *1E10) , H[2][1]* ((3.405E-10) *1E10) , H[2][2]* ((3.405E-10) *1E10) ); printf ("Lagrangian strain norm2 >= NTS_MAX_STRAIN2 " "= %f, exit.\n", (3* (( 0.5 )*( 0.5 )) ) ); save_tmp_config(); exit(1); } } matmul (HIOLD, H, A); matran (A, B); matmul (A, B, G); diag3 (G, x, GA); for (i=0; i<3; i++) { x[i] = (x[i]-1) / 2.; if (fabs(x[i]) >= MAX_STRAIN) { printf ("|strain| in (%.4f %.4f %.4f) direction " "is now greater than\n", GA[i][0], GA[i][1], GA[i][2]); printf ("MAX_STRAIN = %f, repartition the system...\n", MAX_STRAIN); for (j=0; j=kstart) && ((kb-kstart)%kw_gr==0) ) accumulate_and_save_gr ("gr.m" ); for (i=0; i 2) factorHT = 2; } pair_potential(); Tnow = kine * (119.8* 1.380658E-23 ) / neff / (1.5 * 1.380658E-23 ); factorT = 1.-delta/2./tau*((Tnow+1.)/(T+1.)-1); if (factorT < 0.5) factorT = 0.5; else if (factorT > 2) factorT = 2; if (kb >= 1) { kinetum += kine; k1tum += 1./kine; k2tum += 1./kine/kine; for (i=0; i<3; i++) for (j=0; j<3; j++) { Htum[i][j] += H[i][j]; strestum[i][j] += stress[i][j]; v1k1tum[i][j] += virial[i][j] / kine * sqrt(volume); v1k2tum[i][j] += virial[i][j] / kine / kine * sqrt(volume); for (k=0; k<3; k++) for (l=0; l<3; l++) v2k2tum[i][j][k][l] += virial[i][j] * virial[k][l] / kine / kine * volume; } for (i=0; i<6; i++) for (j=0; j<6; j++) ck1tum[i][j] += c[i][j] / kine; } if (kb >= kstart) { potesum += pote; kinesum += kine; k1sum += 1./kine; k2sum += 1./kine/kine; for (i=0; i<3; i++) for (j=0; j<3; j++) { Hsum[i][j] += H[i][j]; stressum[i][j] += stress[i][j]; v1k1sum[i][j] += virial[i][j] / kine * sqrt(volume); v1k2sum[i][j] += virial[i][j] / kine / kine * sqrt(volume); for (k=0; k<3; k++) for (l=0; l<3; l++) v2k2sum[i][j][k][l] += virial[i][j] * virial[k][l] / kine / kine * volume; } for (i=0; i<6; i++) for (j=0; j<6; j++) ck1sum[i][j] += c[i][j] / kine; } if ( (kb>0) && (kb%kw_lp==0) ) { kinetum /= kw_lp; k1tum /= kw_lp; k2tum /= kw_lp; cc = kinetum / (1.5*neff); dd = matinv(Htum,A) / kw_lp / kw_lp / kw_lp; for (i=0; i<6; i++) for (j=0; j<6; j++) ck1tum[i][j] *= (1.5*neff-1) * cc / kw_lp; ck1tum[0][0] += 2*neff*cc/dd; ck1tum[1][1] += 2*neff*cc/dd; ck1tum[2][2] += 2*neff*cc/dd; ck1tum[3][3] += neff*cc/dd; ck1tum[4][4] += neff*cc/dd; ck1tum[5][5] += neff*cc/dd; dE2 = (1.5*neff-1)*(1.5*neff-2)*k2tum - (1.5*neff-1)*(1.5*neff-1)*k1tum*k1tum; for (i=0; i<3; i++) for (j=0; j<3; j++) { A[i][j] = -(1.5*neff-1)*(1.5*neff-1)* k1tum*v1k1tum[i][j]/kw_lp + (1.5*neff-1)*(1.5*neff-2)* v1k2tum[i][j]/kw_lp; v1k1tum[i][j] *= (1.5*neff-1)*sqrt(cc)/kw_lp; for (k=0; k<3; k++) for (l=0; l<3; l++) v2k2tum[i][j][k][l] *= (1.5*neff-1)*(1.5*neff-2)*cc/kw_lp; } ck1tum[0][0] += v1k1tum[0][0] * v1k1tum[0][0] + A[0][0]*A[0][0]*cc/dE2 - v2k2tum[0][0][0][0]; ck1tum[0][1] += v1k1tum[0][0] * v1k1tum[1][1] + A[0][0]*A[1][1]*cc/dE2 - v2k2tum[0][0][1][1]; ck1tum[0][2] += v1k1tum[0][0] * v1k1tum[2][2] + A[0][0]*A[2][2]*cc/dE2 - v2k2tum[0][0][2][2]; ck1tum[1][1] += v1k1tum[1][1] * v1k1tum[1][1] + A[1][1]*A[1][1]*cc/dE2 - v2k2tum[1][1][1][1]; ck1tum[1][2] += v1k1tum[1][1] * v1k1tum[2][2] + A[1][1]*A[2][2]*cc/dE2 - v2k2tum[1][1][2][2]; ck1tum[2][2] += v1k1tum[2][2] * v1k1tum[2][2] + A[2][2]*A[2][2]*cc/dE2 - v2k2tum[2][2][2][2]; ck1tum[0][3] += v1k1tum[0][0] * v1k1tum[1][2] + A[0][0]*A[1][2]*cc/dE2 - v2k2tum[0][0][1][2]; ck1tum[0][4] += v1k1tum[0][0] * v1k1tum[0][2] + A[0][0]*A[0][2]*cc/dE2 - v2k2tum[0][0][0][2]; ck1tum[0][5] += v1k1tum[0][0] * v1k1tum[0][1] + A[0][0]*A[0][1]*cc/dE2 - v2k2tum[0][0][0][1]; ck1tum[1][3] += v1k1tum[1][1] * v1k1tum[1][2] + A[1][1]*A[1][2]*cc/dE2 - v2k2tum[1][1][1][2]; ck1tum[1][4] += v1k1tum[1][1] * v1k1tum[0][2] + A[1][1]*A[0][2]*cc/dE2 - v2k2tum[1][1][0][2]; ck1tum[1][5] += v1k1tum[1][1] * v1k1tum[0][1] + A[1][1]*A[0][1]*cc/dE2 - v2k2tum[1][1][0][1]; ck1tum[2][3] += v1k1tum[2][2] * v1k1tum[1][2] + A[2][2]*A[1][2]*cc/dE2 - v2k2tum[2][2][1][2]; ck1tum[2][4] += v1k1tum[2][2] * v1k1tum[0][2] + A[2][2]*A[0][2]*cc/dE2 - v2k2tum[2][2][0][2]; ck1tum[2][5] += v1k1tum[2][2] * v1k1tum[0][1] + A[2][2]*A[0][1]*cc/dE2 - v2k2tum[2][2][0][1]; ck1tum[3][3] += v1k1tum[1][2] * v1k1tum[1][2] + A[1][2]*A[1][2]*cc/dE2 - v2k2tum[1][2][1][2]; ck1tum[3][4] += v1k1tum[1][2] * v1k1tum[0][2] + A[1][2]*A[0][2]*cc/dE2 - v2k2tum[1][2][0][2]; ck1tum[3][5] += v1k1tum[1][2] * v1k1tum[0][1] + A[1][2]*A[0][1]*cc/dE2 - v2k2tum[1][2][0][1]; ck1tum[4][4] += v1k1tum[0][2] * v1k1tum[0][2] + A[0][2]*A[0][2]*cc/dE2 - v2k2tum[0][2][0][2]; ck1tum[4][5] += v1k1tum[0][2] * v1k1tum[0][1] + A[0][2]*A[0][1]*cc/dE2 - v2k2tum[0][2][0][1]; ck1tum[5][5] += v1k1tum[0][1] * v1k1tum[0][1] + A[0][1]*A[0][1]*cc/dE2 - v2k2tum[0][1][0][1]; for (i=0; i<6; i++) for (j=i+1; j<6; j++) ck1tum[j][i] = ck1tum[i][j]; pressure = (strestum[0][0]+strestum[1][1]+ strestum[2][2])/kw_lp/3.; hamiltonian = pote + kine; if ( scale_H ) hamiltonian += kinewall - hydro*volume - volumeold*(sout[0][0]*(G[0][0]-1)/2 + sout[1][1]*(G[1][1]-1)/2 + sout[2][2]*(G[2][2]-1)/2 + sout[0][1]*G[0][1] + sout[0][2]*G[0][2] + sout[1][2]*G[1][2]); fprintf (property, "%7d %.3f %.6f %.6f %.5f\n", kb, Tnow, pote/np, hamiltonian/np, pressure); fflush (property); if (kb >= kstart) save_mesh ("me.out" ); printf ("\nKB = %d, T_md = %f K, T_real = %f K\n", kb, Tnow, T_md_to_T_real(Tnow)); printf ("Potential energy = %.5f\\%.5f eV" "\\%.4f K\\%.1f J/mol\n", pote/np, pote/np* (119.8* 1.380658E-23 ) / 1.60217733E-19 , pote/np* (119.8* 1.380658E-23 ) / 1.380658E-23 , pote/np* (119.8* 1.380658E-23 ) * 6.0221367E+23 ); printf ("Hamiltonian = %f\\%f eV\\%f K\\%.1f J/mol\n", hamiltonian/np, hamiltonian/np* (119.8* 1.380658E-23 ) / 1.60217733E-19 , hamiltonian/np* (119.8* 1.380658E-23 ) / 1.380658E-23 , hamiltonian/np* (119.8* 1.380658E-23 ) * 6.0221367E+23 ); printf ("System pressure = %f\\%f MPa\n", pressure, pressure* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); printf ("Density = %f\\%f 1/A^3\\%f mol/m^3\n\n", np/volume, np/volume/ pow((3.405E-10) ,3.) *1E-30, np/volume/ 6.0221367E+23 / pow((3.405E-10) ,3.) ); printf (" | %.3f %.3f %.3f |", Htum[0][0]* ((3.405E-10) *1E10) /kw_lp, Htum[0][1]* ((3.405E-10) *1E10) /kw_lp, Htum[0][2]* ((3.405E-10) *1E10) /kw_lp); printf (" | %.3f %.3f %.3f |\n", strestum[0][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) /kw_lp, strestum[0][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) /kw_lp, strestum[0][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) /kw_lp); printf ("H(A) = | %.3f %.3f %.3f |", Htum[1][0]* ((3.405E-10) *1E10) /kw_lp, Htum[1][1]* ((3.405E-10) *1E10) /kw_lp, Htum[1][2]* ((3.405E-10) *1E10) /kw_lp); printf (" t(MPa) = | %.3f %.3f %.3f |\n", strestum[1][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) /kw_lp, strestum[1][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) /kw_lp, strestum[1][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) /kw_lp); printf (" | %.3f %.3f %.3f |", Htum[2][0]* ((3.405E-10) *1E10) /kw_lp, Htum[2][1]* ((3.405E-10) *1E10) /kw_lp, Htum[2][2]* ((3.405E-10) *1E10) /kw_lp); printf (" | %.3f %.3f %.3f |\n\n", strestum[2][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) /kw_lp, strestum[2][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) /kw_lp, strestum[2][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) /kw_lp); printf (" | %8.1f %8.1f %8.1f %8.1f %8.1f %8.1f |\n", ck1tum[0][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[0][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1tum[0][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[0][3]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1tum[0][4]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[0][5]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); printf (" | %8.1f %8.1f %8.1f %8.1f %8.1f %8.1f |\n", ck1tum[1][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[1][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1tum[1][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[1][3]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1tum[1][4]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[1][5]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); printf ("C(MPa) = | %8.1f %8.1f %8.1f %8.1f %8.1f %8.1f |\n", ck1tum[2][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[2][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1tum[2][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[2][3]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1tum[2][4]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[2][5]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); printf (" | %8.1f %8.1f %8.1f %8.1f %8.1f %8.1f |\n", ck1tum[3][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[3][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1tum[3][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[3][3]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1tum[3][4]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[3][5]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); printf (" | %8.1f %8.1f %8.1f %8.1f %8.1f %8.1f |\n", ck1tum[4][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[4][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1tum[4][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[4][3]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1tum[4][4]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[4][5]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); printf (" | %8.1f %8.1f %8.1f %8.1f %8.1f %8.1f |\n\n", ck1tum[5][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[5][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1tum[5][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[5][3]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1tum[5][4]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ,ck1tum[5][5]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); print_list_statistics (((_IO_FILE*)(&_IO_stdout_)) ); printf ("Re-anchoring rate = %f /atom/step.\n\n", (double)anchortum/kw_lp/np); for (cc=dd=i=0; i 0) && (kb % tmp_config_freq == 0) ) save_tmp_config(); for (i=0; i 0) H1[i][j] *= factorHT; } volume = matinv (H, HI); if ( (kb >= Havg_start) && (kb < kstart) ) matadd (H, Havg, Havg); } if ( NTS && (SCALM<=kstart) && (TCALM<=kstart) && (kb==kstart) && (kstart > Havg_start) ) { for (i=0; i<3; i++) for (j=0; j<3; j++) { H[i][j] = Havg[i][j] / (kstart-Havg_start); H1[i][j] = 0.; H2[i][j] = 0.; H3[i][j] = 0.; H4[i][j] = 0.; H5[i][j] = 0.; } volume = matinv (H, HI); } for (i=0; i1.) s[i]--; if ( scale_T && (tau>0) ) s1[i] *= factorT; } if ( CALC && (kb>=kstart) ) { for (i=0; i<3; i++) { cj[i] *= sqrt(((119.8* 1.380658E-23 ) / ((3.405E-10) *sqrt((39.948*1E-3/ 6.0221367E+23 ) / (119.8* 1.380658E-23 ) )) / (3.405E-10) ) *dTdR* delta/(1.380658E-23 *T/ (119.8* 1.380658E-23 ) *T*volume)); fwrite (&cj[i], sizeof(double), 1, heat_current[i]); for (j=0; j<3; j++) stress[i][j] *= sqrt(((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) * ((3.405E-10) *sqrt((39.948*1E-3/ 6.0221367E+23 ) / (119.8* 1.380658E-23 ) )) * delta*volume/(1.380658E-23 *T/ (119.8* 1.380658E-23 ) )); } fwrite (&stress[1][2],sizeof(double),1,shear_stress[0]); fwrite (&stress[0][2],sizeof(double),1,shear_stress[1]); fwrite (&stress[0][1],sizeof(double),1,shear_stress[2]); if ( kb % kw_lp == 0 ) { fflush(heat_current[0]); fflush(heat_current[1]); fflush(heat_current[2]); fflush(shear_stress[0]); fflush(shear_stress[1]); fflush(shear_stress[2]); } } } tmp = (double) (maxkb - kstart); potesum /= tmp; kinesum /= tmp; k1sum /= tmp; k2sum /= tmp; for (i=0; i<3; i++) for (j=0; j<3; j++) { Hsum[i][j] /= tmp; stressum[i][j] /= tmp; v1k1sum[i][j] /= tmp; v1k2sum[i][j] /= tmp; for (k=0; k<3; k++) for (l=0; l<3; l++) v2k2sum[i][j][k][l] /= tmp; } for (i=0; i<6; i++) for (j=0; j<6; j++) ck1sum[i][j] /= tmp; volume = matinv (Hsum, HI); cc = kinesum / (1.5*neff); Tnow = cc * (119.8* 1.380658E-23 ) / 1.380658E-23 ; Treal = T_md_to_T_real (Tnow); T_real_to_T_md (Treal, &dTdR); dE2 = (1.5*neff-1)*(1.5*neff-2)*k2sum - (1.5*neff-1)*(1.5*neff-1)*k1sum*k1sum; dV2 = dEdV = dd = tmp = 0.; for (i=0; i<3; i++) { dV2 += v1k1sum[i][i]/3./sqrt(volume); dEdV += (1.5*neff-1)*(1.5*neff-2)*v1k2sum[i][i]/3./sqrt(volume) - (1.5*neff-1)*(1.5*neff-1)*v1k1sum[i][i]*k1sum/3./sqrt(volume); for (k=0; k<3; k++) { dd += v2k2sum[i][i][k][k]; tmp += ck1sum[i][k]; } } dV2 = -neff/volume/volume -(1.5*neff-1)*(1.5*neff-1)*dV2*dV2 +(1.5*neff-1)*(1.5*neff-2)*dd/9./volume -2.*(1.5*neff-1)/3./volume*dV2 -(1.5*neff-1)/9./volume*(tmp-3*dV2); pressure = (stressum[0][0]+stressum[1][1]+stressum[2][2])/3.; bulk_modulus = cc*volume*(dEdV*dEdV/dE2-dV2); alpha = dTdR*(pressure-dEdV/dE2)/Tnow/bulk_modulus/3.; Cv = -(1.380658E-23 *T/ (119.8* 1.380658E-23 ) )/T*(1.5*neff-1)*(1.5*neff-1)* k1sum*k1sum/dE2*dTdR; Cp = Cv+9.*volume*Treal*bulk_modulus*alpha*alpha; for (i=0; i<6; i++) for (j=0; j<6; j++) ck1sum[i][j] *= (1.5*neff-1)*cc; ck1sum[0][0] += 2*neff*cc/volume; ck1sum[1][1] += 2*neff*cc/volume; ck1sum[2][2] += 2*neff*cc/volume; ck1sum[3][3] += neff*cc/volume; ck1sum[4][4] += neff*cc/volume; ck1sum[5][5] += neff*cc/volume; for (i=0;i<3;i++) for (j=0;j<3;j++) { A[i][j] = -(1.5*neff-1)*(1.5*neff-1)*k1sum*v1k1sum[i][j] +(1.5*neff-1)*(1.5*neff-2)*v1k2sum[i][j]; v1k1sum[i][j] *= (1.5*neff-1)*sqrt(cc); for (k=0;k<3;k++) for (l=0;l<3;l++) v2k2sum[i][j][k][l] *= (1.5*neff-1)*(1.5*neff-2)*cc; } ck1sum[0][0] += v1k1sum[0][0]*v1k1sum[0][0]+A[0][0]*A[0][0]*cc/dE2 -v2k2sum[0][0][0][0]; ck1sum[0][1] += v1k1sum[0][0]*v1k1sum[1][1]+A[0][0]*A[1][1]*cc/dE2 -v2k2sum[0][0][1][1]; ck1sum[0][2] += v1k1sum[0][0]*v1k1sum[2][2]+A[0][0]*A[2][2]*cc/dE2 -v2k2sum[0][0][2][2]; ck1sum[1][1] += v1k1sum[1][1]*v1k1sum[1][1]+A[1][1]*A[1][1]*cc/dE2 -v2k2sum[1][1][1][1]; ck1sum[1][2] += v1k1sum[1][1]*v1k1sum[2][2]+A[1][1]*A[2][2]*cc/dE2 -v2k2sum[1][1][2][2]; ck1sum[2][2] += v1k1sum[2][2]*v1k1sum[2][2]+A[2][2]*A[2][2]*cc/dE2 -v2k2sum[2][2][2][2]; ck1sum[0][3] += v1k1sum[0][0]*v1k1sum[1][2]+A[0][0]*A[1][2]*cc/dE2 -v2k2sum[0][0][1][2]; ck1sum[0][4] += v1k1sum[0][0]*v1k1sum[0][2]+A[0][0]*A[0][2]*cc/dE2 -v2k2sum[0][0][0][2]; ck1sum[0][5] += v1k1sum[0][0]*v1k1sum[0][1]+A[0][0]*A[0][1]*cc/dE2 -v2k2sum[0][0][0][1]; ck1sum[1][3] += v1k1sum[1][1]*v1k1sum[1][2]+A[1][1]*A[1][2]*cc/dE2 -v2k2sum[1][1][1][2]; ck1sum[1][4] += v1k1sum[1][1]*v1k1sum[0][2]+A[1][1]*A[0][2]*cc/dE2 -v2k2sum[1][1][0][2]; ck1sum[1][5] += v1k1sum[1][1]*v1k1sum[0][1]+A[1][1]*A[0][1]*cc/dE2 -v2k2sum[1][1][0][1]; ck1sum[2][3] += v1k1sum[2][2]*v1k1sum[1][2]+A[2][2]*A[1][2]*cc/dE2 -v2k2sum[2][2][1][2]; ck1sum[2][4] += v1k1sum[2][2]*v1k1sum[0][2]+A[2][2]*A[0][2]*cc/dE2 -v2k2sum[2][2][0][2]; ck1sum[2][5] += v1k1sum[2][2]*v1k1sum[0][1]+A[2][2]*A[0][1]*cc/dE2 -v2k2sum[2][2][0][1]; ck1sum[3][3] += v1k1sum[1][2]*v1k1sum[1][2]+A[1][2]*A[1][2]*cc/dE2 -v2k2sum[1][2][1][2]; ck1sum[3][4] += v1k1sum[1][2]*v1k1sum[0][2]+A[1][2]*A[0][2]*cc/dE2 -v2k2sum[1][2][0][2]; ck1sum[3][5] += v1k1sum[1][2]*v1k1sum[0][1]+A[1][2]*A[0][1]*cc/dE2 -v2k2sum[1][2][0][1]; ck1sum[4][4] += v1k1sum[0][2]*v1k1sum[0][2]+A[0][2]*A[0][2]*cc/dE2 -v2k2sum[0][2][0][2]; ck1sum[4][5] += v1k1sum[0][2]*v1k1sum[0][1]+A[0][2]*A[0][1]*cc/dE2 -v2k2sum[0][2][0][1]; ck1sum[5][5] += v1k1sum[0][1]*v1k1sum[0][1]+A[0][1]*A[0][1]*cc/dE2 -v2k2sum[0][1][0][1]; for (i=0; i<6; i++) for (j=i+1; j<6; j++) ck1sum[j][i] = ck1sum[i][j]; fprintf (property, "\n This is an MD simulation of " "%d particles using\n", np); fprintf (property, " %s pair potential with RCUT = %.3f " "= %.4f A.\n\n", potfun_name, (2.5* (3.405E-10) / (3.405E-10) ) , (2.5* (3.405E-10) / (3.405E-10) ) * ((3.405E-10) *1E10) ); if (nkind>1) print_atom_statistics(property, " "); fprintf (property, " During the entire simulation (%d x %.4f ps),\n", maxkb, delta* (((3.405E-10) *sqrt((39.948*1E-3/ 6.0221367E+23 ) / (119.8* 1.380658E-23 ) )) *1E12) ); fprintf (property, " there are %d strain-triggered repartitions" " (MAX_STRAIN = %.3f),\n", nw_repartition, MAX_STRAIN); fprintf (property, " and %d re-anchors, or %f /atom/step.\n\n", nw_anchor, (double)nw_anchor/np/maxkb); fprintf (property, " The final average starts from step " "%d to %d.\n\n", kstart, maxkb); fprintf (property, " The average MD temperature is %f = %f K,\n", Tnow/T*(1.380658E-23 *T/ (119.8* 1.380658E-23 ) ), Tnow); fprintf (property, " corresponding to a real temperature of %f K,\n", Treal); fprintf (property, " with dT_md / dT_real = %f.\n\n", dTdR); fprintf (property, " H[1][1] = %8.4f A stress[1][1] = %9.4f MPa\n", Hsum[0][0]* ((3.405E-10) *1E10) , stressum[0][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); fprintf (property, " H[1][2] = %8.4f A stress[1][2] = %9.4f MPa\n", Hsum[0][1]* ((3.405E-10) *1E10) , stressum[0][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); fprintf (property, " H[1][3] = %8.4f A stress[1][3] = %9.4f MPa\n", Hsum[0][2]* ((3.405E-10) *1E10) , stressum[0][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); fprintf (property, " H[2][1] = %8.4f A stress[2][1] = %9.4f MPa\n", Hsum[1][0]* ((3.405E-10) *1E10) , stressum[1][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); fprintf (property, " H[2][2] = %8.4f A stress[2][2] = %9.4f MPa\n", Hsum[1][1]* ((3.405E-10) *1E10) , stressum[1][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); fprintf (property, " H[2][3] = %8.4f A stress[2][3] = %9.4f MPa\n", Hsum[1][2]* ((3.405E-10) *1E10) , stressum[1][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); fprintf (property, " H[3][1] = %8.4f A stress[3][1] = %9.4f MPa\n", Hsum[2][0]* ((3.405E-10) *1E10) , stressum[2][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); fprintf (property, " H[3][2] = %8.4f A stress[3][2] = %9.4f MPa\n", Hsum[2][1]* ((3.405E-10) *1E10) , stressum[2][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); fprintf (property, " H[3][3] = %8.4f A stress[3][3] = %9.4f MPa\n\n", Hsum[2][2]* ((3.405E-10) *1E10) , stressum[2][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); fprintf (property, " Avg. density = %f\\%f 1/A^3\\%f mol/m^3\n", np/volume, np/volume/ pow((3.405E-10) ,3.) *1E-30, np/volume/ 6.0221367E+23 / pow((3.405E-10) ,3.) ); fprintf (property, " Avg. pressure = %f\\%f MPa\n\n", pressure, pressure* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); fprintf (property, " Avg. kinetic energy = " "%.4f\\%.4f eV\\%.1f K\\%.1f J/mol\n", kinesum/np, kinesum/np* (119.8* 1.380658E-23 ) / 1.60217733E-19 , kinesum/np* (119.8* 1.380658E-23 ) / 1.380658E-23 , kinesum* (119.8* 1.380658E-23 ) * 6.0221367E+23 /np); fprintf (property, " Avg. potential energy = " "%.4f\\%.4f eV\\%.1f K\\%.1f J/mol\n", potesum/np, potesum/np* (119.8* 1.380658E-23 ) / 1.60217733E-19 , potesum/np* (119.8* 1.380658E-23 ) / 1.380658E-23 , potesum* (119.8* 1.380658E-23 ) * 6.0221367E+23 /np); potesum += kinesum; fprintf (property, " Avg. internal energy = " "%.4f\\%.4f eV\\%.1f K\\%.1f J/mol\n", potesum/np, potesum/np* (119.8* 1.380658E-23 ) / 1.60217733E-19 , potesum/np* (119.8* 1.380658E-23 ) / 1.380658E-23 , potesum* (119.8* 1.380658E-23 ) * 6.0221367E+23 /np); potesum += pressure * volume; fprintf (property, " Avg. enthalpy = " "%.4f\\%.4f eV\\%.1f K\\%.1f J/mol\n\n", potesum/np, potesum/np* (119.8* 1.380658E-23 ) / 1.60217733E-19 , potesum/np* (119.8* 1.380658E-23 ) / 1.380658E-23 , potesum* (119.8* 1.380658E-23 ) * 6.0221367E+23 /np); fprintf (property, " | %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f |\n", ck1sum[0][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[0][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[0][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[0][3]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[0][4]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[0][5]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); fprintf (property, " | %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f |\n", ck1sum[1][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[1][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[1][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[1][3]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[1][4]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[1][5]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); fprintf (property, " C(MPa) = | %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f |\n", ck1sum[2][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[2][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[2][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[2][3]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[2][4]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[2][5]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); fprintf (property, " | %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f |\n", ck1sum[3][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[3][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[3][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[3][3]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[3][4]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[3][5]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); fprintf (property, " | %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f |\n", ck1sum[4][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[4][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[4][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[4][3]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[4][4]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[4][5]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); fprintf (property, " | %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f |\n\n", ck1sum[5][0]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[5][1]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[5][2]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[5][3]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[5][4]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) , ck1sum[5][5]* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); fprintf (property, " The finite stress bulk modulus " "BT = (c11+2*c12+p)/3 = %f MPa\n", bulk_modulus* (((119.8* 1.380658E-23 ) / pow((3.405E-10) ,3.) ) /1.E6) ); fprintf (property, " The line thermal expansion coefficient alpha = %f/K\n", alpha); fprintf (property, " The heat capacity Cv = %f J/g/K = %f J/mol/K\n", Cv* (119.8* 1.380658E-23 ) /(totalmass* (39.948*1E-3/ 6.0221367E+23 ) *1000.), Cv* (119.8* 1.380658E-23 ) /neff* 6.0221367E+23 ); fprintf (property, " Cp = %f J/g/K = %f J/mol/K\n", Cp* (119.8* 1.380658E-23 ) /(totalmass* (39.948*1E-3/ 6.0221367E+23 ) *1000.), Cp* (119.8* 1.380658E-23 ) /neff* 6.0221367E+23 ); fprintf (property, "\n The self-diffusion coefficient D = %f 10^(-10) m^2/s\n", (msd-msd_kstart)* (3.405E-10) * (3.405E-10) / ((maxkb-kstart-1)*delta* ((3.405E-10) *sqrt((39.948*1E-3/ 6.0221367E+23 ) / (119.8* 1.380658E-23 ) )) )/6.*1E10); fclose (structure); write_config (fn_config_write); free (idx); free (list); free (bindx); free (binlist); free (binbindx); free (binbinlist); free (mb); free (tpd); free (tp[0]); free (tp); free (mybin); free (ep); free (mass); free (d); free (fs); free (f); free (v); free (s5); free (s4); free (s3); free (s2); free (s1); free (s); cancel_all_schedules(); if ( CALC ) { for (i=0; i<3; i++) { fclose (heat_current[i]); fclose (shear_stress[i]); } smile (0 , property); smile (1 , property); } cc = watch ("mission", 3); fprintf (property, "\n This simulation took %s --" "\n or %f s/step, %f ms/step/particle.\n\n", earth_time(cc), cc/maxkb, 1E3*cc/np/maxkb); fclose (property); return(0); }