Descrever as coordenadas de um sistemas requer definir um centro, ou uma origem. A escolha dessa origem pode ser feita por um questão gosto do usuário ou conveniência de cálculos. Na Física, somos aconselhados a usar o Centro de Massas (CM) como origem do nosso sistemas de coordenadas pois qualquer distribuição de massa (seja um conjunto de partículas como uma molécula ou um sistema maciço) pode ser tratado como se toda a massa estivesse localizada no CM. Isso facilita no estudo da cinemática/dinâmica de todo o sistema, mesmo que internamente não faça diferença a posição da origem. Trabalhando no CM garantimos que o momento angular do sistema se conserve. Em simulações é importante para aplicar correções de pressão e temperatura.
Entretanto não é trivial ter um programa que faça o cálculo do CM e o deslocamento para o mesmo. Há algumas ferramentas já fechadas na qual não é possível alterar o código a gosto do usuário. A proposta aqui é fornecer um código acessível aos usuários que trabalham com moléculas, que possa ser usado do jeito que foi fornecido ou modificado.
Segue o programa que escrevi, em Fortran 90 para calcular o centro de massas de um arquivo XYZ. Ele também desloca a origem do sistema de coordenadas do arquivo de entrada para o centro de massas. O código é livre para fazer qualquer alteração. Note que há algumas opções de argumentos de comandos e alocação de memória, que são interessanes para quem está estudando o Fortran.
PROGRAM xyzcom
! =====================================================================
! Program xyzcom
! - calculate the center of mass (COM) from a xyz file
! - change the origin on initial XYZ file to COM
!
! Lucas Modesto da Costa - lucmod@gmail.com
! 20/01/2020 -
IMPLICIT NONE
CHARACTER(len=40) :: filein,fileout, title
CHARACTER(len=40),ALLOCATABLE,DIMENSION(:) :: args
CHARACTER(len=2) :: findatom
INTEGER :: natom,ngroup,ni,nj
INTEGER :: i,j
INTEGER :: ix,num_args
REAL, ALLOCATABLE, DIMENSION(:) :: x, y, z
REAL :: xcm=0, ycm=0.0, zcm=0, tmass=0
INTEGER, ALLOCATABLE, DIMENSION(:) :: indx
CHARACTER(len=5),ALLOCATABLE,DIMENSION(:) :: name
REAL, DIMENSION(105) :: atommass
CHARACTER(len=2),DIMENSION(105) :: atomname
LOGICAL :: located, fileexist
!List of atom mass
DATA atommass &
/ 1.008, 4.003, 6.939, 9.012, 10.811, 12.011, 14.007, &
15.999, 18.998, 20.183, 22.989, 24.312, 26.982, 28.086, &
30.974, 32.064, 35.453, 39.948, 39.102, 40.080, 44.956, &
47.900, 50.942, 51.996, 54.938, 55.847, 58.933, 58.710, &
63.540, 65.370, 69.720, 72.590, 74.922, 78.960, 79.909, &
83.800, 85.470, 87.620, 88.905, 91.220, 92.906, 95.940, &
98.000, 101.070, 102.905, 106.400, 107.870, 112.400, 114.820, &
118.690, 121.750, 127.600, 126.904, 131.300, 132.905, 137.340, &
138.910, 140.120, 140.907, 144.240, 147.000, 150.350, 151.960, &
157.250, 158.924, 162.500, 164.930, 167.260, 168.934, 173.040, &
174.970, 178.490, 180.948, 183.850, 186.200, 190.200, 192.200, &
195.090, 196.967, 200.590, 204.370, 207.190, 208.980, 210.000, &
210.000, 222.000, 223.000, 226.000, 227.000, 232.038, 231.000, &
238.030, 237.000, 242.000, 243.000, 247.000, 247.000, 249.000, &
254.000, 253.000, 256.000, 254.000, 257.000, 0.000, 0.000 /
!List of atom name
DATA atomname &
/ ' H', 'He', 'Li', 'Be', ' B', ' C', ' N', ' O', ' F', 'Ne', 'Na',&
'Mg', 'Al', 'Si', ' P', ' S', 'Cl', 'Ar', ' K', 'Ca', 'Sc', 'Ti',&
' V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As',&
'Se', 'Br', 'Kr', 'Rb', 'Sr', ' Y', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru',&
'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', ' I', 'Xe', 'Cs',&
'Ba', 'La', 'Ce', 'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy',&
'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', ' W', 'Re', 'Os', 'Ir',&
'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra',&
'Ac', 'Th', 'Pa', ' U', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es',&
'Fm', 'Md', 'No', 'Lw', 'Xx', 'Bq'/
! number of arguments in terminal command
num_args = command_argument_count()
IF (num_args.eq.0) THEN
write(*,*) "No argument has found. "
write(*,*) "Usage: xyzcom <input_file.xyz> [<output_file.xyz>]"
STOP
ENDIF
ALLOCATE(args(num_args))
IF (num_args.gt.2) THEN
write(*,*) "Several arguments has found. "
write(*,*) "Usage: xyzcom <input_file.xyz> [<output_file.xyz>]"
STOP
ENDIF
! read the arguments
DO ix = 1, num_args
CALL get_command_argument(ix,args(ix))
! now parse the argument as you wish
END DO
filein = TRIM(args(1))
INQUIRE( FILE=filein, EXIST=fileexist )
IF ( fileexist) THEN
write(*,*) "Reading ", TRIM(filein)
ELSE
write(*,*) "File ", TRIM(filein), " does not exist."
STOP
ENDIF
OPEN( UNIT= 3, FILE= filein, STATUS="old" )
READ(3,*) natom
ALLOCATE( x(natom), y(natom), z(natom), indx(natom), name(natom))! alocacao de memoria para bw
READ(3,*) title
DO i=1, natom
READ(3,*) name(i),x(i),y(i),z(i)
WRITE(findatom,'(2A)') TRIM(name(i))
IF (LEN_TRIM(name(i)).eq.1) THEN
WRITE(findatom,'(2A)') " ",TRIM(name(i))
ELSEIF (LEN_TRIM(name(i)).ge.3) THEN
WRITE(*,*) "ERROR: ", name(i), " not supported. Maybe the 3 or more letters to a atom name is a mistake."
STOP
ELSE
ENDIF
IF (ANY(INDEX(atomname, findatom).eq.1)) THEN
indx(i) = MAXLOC(INDEX(atomname, findatom,located),1)
ELSE
WRITE(*,*) "ERROR: atom ",findatom," is not supported. Please pay attetion correct name or format of then atom."
STOP
ENDIF
ENDDO
CLOSE( unit= 3 )
DO i=1, natom
! WRITE(*,'(A2,2X,F6.3,3(2X,F9.6))') atomname(indx(i)), atommass(indx(i)),x(i),y(i),z(i)
tmass = tmass+ atommass(indx(i))
xcm = xcm + x(i)*atommass(indx(i))
ycm = ycm + y(i)*atommass(indx(i))
zcm = zcm + z(i)*atommass(indx(i))
ENDDO
xcm = xcm/tmass
ycm = ycm/tmass
zcm = zcm/tmass
WRITE(*,*)"Total mass - Coordinates center of mass"
WRITE(*,'(4(2X,F15.9))')tmass,xcm, ycm, zcm
!
! DO i=1, natom
! WRITE(*,'(A2,2X,F6.3,3(2X,F9.6))') atomname(indx(i)), atommass(indx(i)),x(i)-xcm,y(i)-ycm,z(i)-zcm
! ENDDO
IF (num_args.eq.2) THEN
fileout = TRIM(args(2))
OPEN( UNIT= 4, FILE= fileout, STATUS="UNKNOWN" )
WRITE(*,*) "Write coordinates with center of mass origin in ", TRIM(fileout)
WRITE(4,*) natom
WRITE(4,*) title
DO i=1, natom
WRITE(4,'(A2,4X,3(F15.9,2X))') atomname(indx(i)), x(i)-xcm,y(i)-ycm,z(i)-zcm
ENDDO
ENDIF
ENDPROGRAM xyzcom
Sei que o Fortran tem sido colocado de lado outras linguagens bem mais práticas. Esse código foi escrito e validado em 8 horas, sendo que havia muito tempo que não usava o Fortran. Fiz uma pesquisa rápida na net para uma alternativa em Python. O programa que desempenha a mesma função em python foi criado e testado em 5 minutos. Muito tentador. A fonte desse código está aqui.
# -*- coding: utf-8 -*-
#!/usr/bin/env python
from molmod import *
# 0) Load the molecule from file
mol1 = Molecule.from_file("lixo.xyz")
# 1) The atomic masses are not included in the sdf file. We need them below to
# compute the center of mass, so we assign standard masses.
mol1.set_default_masses()
print(mol1.masses/amu)
# 2) Create modified coordinates. The attribute coordinates is an N by 3 array,
# while the com attribute is an array with three components. The following
# subtraction is performed elementwise at the level of the rows. Each row
# of the coordinates array corresponds to one atom. From each atomic coordinate
# the center of mass is subtracted to construct a new coordinates array.
new_coordinates = mol1.coordinates - mol1.com
# 3) Make a copy of mol1 with updated coordinates and write it to a file.
mol2 = mol1.copy_with(coordinates=new_coordinates)
mol2.write_to_file("lixo_com.xyz")
print("Written file lixo_com.xyz")
Prof. Lucas Modesto da Costa 