Amber-MD / cpptraj

Biomolecular simulation trajectory/data analysis.
Other
132 stars 63 forks source link

Error: No reference atoms selected for parm , [@CA] #1083

Closed H-EKE closed 1 month ago

H-EKE commented 3 months ago

Hi,

I am trying to replicate the analysis that I did with MDAnaalysis on cpptraj

import MDAnalysis as mda
from MDAnalysis.analysis import rms, align
# import nglview as nv

import warnings
# suppress some MDAnalysis warnings about writing PDB files
warnings.filterwarnings('ignore')

u = mda.Universe('b.psf', 'cent.xtc')

average = align.AverageStructure(u, u, select='protein and name CA',
                                 ref_frame=0).run()
ref = average.results.universe

aligner = align.AlignTraj(u, ref,
                          select='protein and name CA',
                          in_memory=True).run()

c_alphas = u.select_atoms('protein and name CA')
R = rms.RMSF(c_alphas).run()

# Write RMSF values to a file
with open('rmsf_values.txt', 'w') as f:
    for residue, r_value in zip(c_alphas.residues, R.results.rmsf):
        f.write(f"{residue} {r_value}\n")

RMSF.in

for i=1;i<21;i++
parm prod1/structure.psf
trajin prod$i/output.xtc
strip :SOD,CLA,TIP3
autoimage 
rms first @CA
average crdset MyAvg$i @CA
rms ref MyAvg$i @CA 
atomicfluct out rmsf_u2.dat @CA byres
go
clear all
done

when I execute RMSF.in to calculate the rmsf , I get the following output

  1: [autoimage]
    Using first molecule as anchor.
    1 molecules are fixed to anchor: 2
    0 molecules are mobile.
  2: [rms first @CA]
    Target mask: [@CA](210)
    Reference topology: structure.psf
    Reference mask: [@CA](210)
  3: [average crdset MyAvg1 @CA]
    Mask [@CA] corresponds to 210 atoms.
    Averaging over 210 atoms.
Warning: Atom selection < total # atoms, stripping parm for averaging only:
  4: [rms ref MyAvg1 @CA]
    Target mask: [@CA](210)
    Reference topology: 
    Reference mask: [@CA](0)
Error: No reference atoms selected for parm , [@CA]
Warning: Setup incomplete for [rms]: Skipping
  5: [atomicfluct out rmsf_u2.dat @CA byres]
    Mask [@CA] corresponds to 210 atoms.

How I can get rid off from

Error: No reference atoms selected for parm , [@CA] 
Warning: Setup incomplete for [rms]: Skipping

Thanks in advance!

drroe commented 3 months ago

Hi. There are two issues. The first is that clear all really will clear everything (including data) so that's overkill. You can just read in the topology once and clear the input trajectory at each loop iteration. The second issue is that you are trying to create an average structure with average, and then use that average structure in the same trajectory pass. You need to first create the average, the align to the average in a second pass. So try something like:

parm prod1/structure.psf
for i=1;i<21;i++
  trajin prod$i/output.xtc
  strip :SOD,CLA,TIP3
  autoimage 
  rms first @CA
  average crdset MyAvg$i @CA
  go
  rms ref MyAvg$i @CA 
  atomicfluct F$i out rmsf_u2.dat @CA byres
  dataset F$i legend Fluct$i
  go
  clear trajin
done

I also added a dataset legend (via dataset legend) to your atomic fluctuation sets to make them easier to distinguish in the output file. Let me know if this works for you.

H-EKE commented 3 months ago

Dear @drroe

Thanks for your fast and helpful reply!

Now i don't have there error but the warning still

Warning: Atom selection < total # atoms, stripping parm for averaging only:
----- output.xtc (1-20000, 1) -----
 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% Complete.

and also the new values for RMSF are extremely high , almost 10 higher that the one repported by MDAnalysis. Why this could happen?

drroe commented 3 months ago

Hi,

I think MDanalysis might report distances as nm, whereas cpptraj uses Angstroms. Maybe this is the source of the 10x discrepancy.

The warning can be safely ignored. It means the topology for your average structure is stripped (i.e. has atoms removed) compared to the input topology, as opposed to the 'strip' command which modifies the topology for all following actions. Hopefully this makes sense. Let me know if not.

H-EKE commented 3 months ago

Thanks for your reply @drroe

MDAnalysis reports distance as Angstroms (https://userguide.mdanalysis.org/stable/units.html).

I think error comes from how I set up the analysis. I wanted just to do a symple rmsf analysis,

so my idea was i 1) fit trajectory to first frame using @CA 2)create average structure using @CA 3) align to that average structure 4 ) And the calculate how fluctuates based on @CA res

drroe commented 3 months ago

Well, the script I posted above should be doing these things. I'm not sure why you're getting different number from MDAnalysis; I'm not very familiar with MDAnalysis syntax, but based on what input you posted, it does look like you are looking at two different trajectories (cent.xtc for MDAnalysis and output.xtc for cpptraj).

The brief test I did with script I posted seemed to perform as expected to give reasonable RMSF values for a stable trajectory. Maybe if you post or attach the entire output from cpptraj that will shed some light on what may be different. Also, if you want to attach your topology and a few trajectory frames I could try to reproduce what you are seeing.

H-EKE commented 3 months ago

output.xtc and cent.xtc are the same trajectories. The difference is thatcent.xtc has already ions, water tripped and is it aligned, to speed up analysis on MDAnalysis. output.log

I have attached the output here.

I tried also with vmd, and I get similar values as MDAnalysis.

File cent.xtc comes from this script

#cent.in
parm prod1/structure.psf
for i=1;i<21;i++
trajin prod$i/output.xtc
strip :SOD,CLA,TIP3
autoimage
rms first @CA 
trajout prod$i/cent.xtc
go
clear trajin
run
done

This script gives the expected output

parm prod1/b.psf
for i=1;i<2;i++
trajin prod$i/cent.xtc 1 last
average crdset MyAvg @CA
run
rms ref MyAvg @CA
atomicfluct out rmsf_u2.dat @CA byres
go
clear trajin
done

My question would be, why if i strip and save the trajectory as in cent.in or if I do it at once as in the initial script, I dont get the same values? Did i miss something?

drroe commented 3 months ago

Sorry for the delay.

My question would be, why if i strip and save the trajectory as in cent.in or if I do it at once as in the initial script, I dont get the same values?

So you're saying that the first script with strip gave problematic values:

parm prod1/structure.psf
for i=1;i<21;i++
  trajin prod$i/output.xtc
  strip :SOD,CLA,TIP3
  autoimage 
  rms first @CA
  average crdset MyAvg$i @CA
  go
  rms ref MyAvg$i @CA 
  atomicfluct F$i out rmsf_u2.dat @CA byres
  dataset F$i legend Fluct$i
  go
  clear trajin
done

But the second set of scripts in which you first strip and image the trajectories and then in a second step calculate atomic fluctuations gave you values that you expect?

drroe commented 3 months ago

I do see one mistake in your output.log. You have an extra run after the rms command but before the average command (I emphasized it with asterisks below):

  BLOCK  0: for ((i=1; i<2; i+=1)) do
  [trajin prod$i/output.xtc 1 last]
    Added command 'trajin' to control block 0.
  [strip :SOD,CLA,TIP3]
    Added command 'strip' to control block 0.
  [autoimage]
    Added command 'autoimage' to control block 0.
  [rms first @CA]
    Added command 'rms' to control block 0.
  **[run]**
    Added command 'run' to control block 0.
  [average crdset MyAvg$i @CA]
    Added command 'average' to control block 0.
  [run]
    Added command 'run' to control block 0.

So what you are doing here is imaging and rms-fitting, then in a second separate pass calculating the average structure without imaging or fitting, i.e. your average structure is not imaged or fit. If you remove the run with asterisks around it I think things should work.

H-EKE commented 3 months ago

That worked perfectly @drroe, thanks a lot again!! I thought that imaging would apply to further analysis, but apparently not. Again thanks a lot for the help!

drroe commented 3 months ago

Glad it worked!

I thought that imaging would apply to further analysis, but apparently not.

So the way that actions work in cpptraj, when you enter one (e.g. rms) it is active for the next run only. Once you type run or go, all currently queued actions (i.e. the stuff printed after ACTION SETUP FOR PARM in the output) are run on all currently loaded trajectories. When the run is finished, those actions are considered completed and are cleared. Trajectories are not automatically cleared to easily allow multiple passes.

If you want to modify coordinates for subsequent passes you can save the coordinates in memory with createcrd, then use crdaction to run additional actions on them, e.g.

parm myparm.parm7
trajin mytraj.nc
rms first @CA 
createcrd MyCrd
run
crdaction MyCrd average @CA crdset MyAvg
crdout MyAvg average.rst7

Hope this helps!

H-EKE commented 3 months ago

It helps a lot! Let me know if there is a way that I can give back the help. Thanks again @drroe

H-EKE commented 1 month ago

Just a small doubt.

I am trying to use a similar script with the mpi version,

parm prod1/structure.psf
for i=1;i<6;i++
  trajin prod$i/output.xtc 1 5000 50
  strip :SOD,CLA,TIP3,POPC
  autoimage 
  rms first @CA
done
average crdset MyAvg @CA
go
rms ref MyAvg @CA
atomicfluct out rmsfA.dat ^3@CA byres
atomicfluct out rmsf465A.dat ^1,@CA byres
run

but I always get this error

CPPTRAJ: Trajectory Analysis. V6.4.4 (GitHub) MPI OpenMP
    ___  ___  ___  ___
     | \/ | \/ | \/ | 
    _|_/\_|_/\_|_/\_|_

| Running on 16 processes.
| 1 OpenMP threads available.
| Date/time: 07/20/24 17:21:37
| Available memory: 262.262 GB

INPUT: Reading input from 'STDIN'
  [parm prod1/structure.psf]
        Reading 'prod1/structure.psf' as Charmm PSF
    Reading Charmm PSF file structure.psf as topology file.
        PSF contains 570226 atoms, 145437 residues.
Warning: Determining bond length parameters from element types for 'structure.psf'.
  [for i=1;i<6;i++]
    Setting up 'for' loop.
CONTROL: Parsing control block.
  BLOCK  0: for ((i=1; i<6; i+=1)) do
  [trajin prod$i/output.xtc 1 5000 50]
        Added command 'trajin' to control block 0.
  [strip :SOD,CLA,TIP3,POPC]
        Added command 'strip' to control block 0.
  [autoimage]
        Added command 'autoimage' to control block 0.
  [rms first @CA]
        Added command 'rms' to control block 0.
  [done]
  BLOCK  0: END
CONTROL: Executing 1 control block(s).
        Loop over 'i' will execute for 5 iterations.
  [trajin prod1/output.xtc 1 5000 50]
        Reading 'prod1/output.xtc' as Gromacs XTC
Warning: Trajectory box type is 'Orthorhombic' but topology box type is 'None'.
Warning: Setting topology box information from trajectory.
  [strip :SOD,CLA,TIP3,POPC]
    STRIP: Stripping atoms in mask [:SOD,CLA,TIP3,POPC]
  [autoimage]
    AUTOIMAGE: To box center based on center of mass, anchor is first molecule.
  [rms first @CA]
    RMSD: (@CA), reference is first frame (@CA).
        Best-fit RMSD will be calculated, coords will be rotated and translated.
  [trajin prod2/output.xtc 1 5000 50]
        Reading 'prod2/output.xtc' as Gromacs XTC
  [strip :SOD,CLA,TIP3,POPC]
    STRIP: Stripping atoms in mask [:SOD,CLA,TIP3,POPC]
  [autoimage]
    AUTOIMAGE: To box center based on center of mass, anchor is first molecule.
  [rms first @CA]
    RMSD: (@CA), reference is first frame (@CA).
        Best-fit RMSD will be calculated, coords will be rotated and translated.
  [trajin prod3/output.xtc 1 5000 50]
        Reading 'prod3/output.xtc' as Gromacs XTC
  [strip :SOD,CLA,TIP3,POPC]
    STRIP: Stripping atoms in mask [:SOD,CLA,TIP3,POPC]
  [autoimage]
    AUTOIMAGE: To box center based on center of mass, anchor is first molecule.
  [rms first @CA]
    RMSD: (@CA), reference is first frame (@CA).
        Best-fit RMSD will be calculated, coords will be rotated and translated.
  [trajin prod4/output.xtc 1 5000 50]
        Reading 'prod4/output.xtc' as Gromacs XTC
  [strip :SOD,CLA,TIP3,POPC]
    STRIP: Stripping atoms in mask [:SOD,CLA,TIP3,POPC]
  [autoimage]
    AUTOIMAGE: To box center based on center of mass, anchor is first molecule.
  [rms first @CA]
    RMSD: (@CA), reference is first frame (@CA).
        Best-fit RMSD will be calculated, coords will be rotated and translated.
  [trajin prod5/output.xtc 1 5000 50]
        Reading 'prod5/output.xtc' as Gromacs XTC
  [strip :SOD,CLA,TIP3,POPC]
    STRIP: Stripping atoms in mask [:SOD,CLA,TIP3,POPC]
  [autoimage]
    AUTOIMAGE: To box center based on center of mass, anchor is first molecule.
  [rms first @CA]
    RMSD: (@CA), reference is first frame (@CA).
        Best-fit RMSD will be calculated, coords will be rotated and translated.
CONTROL: Control block finished.

  [average crdset MyAvg]
        Setting active reference for distance-based masks: 'MyAvg'
    AVERAGE: Averaging over coordinates in mask [*]
        Start: 1  Stop: Final frame
        Saving averaged coords to set 'MyAvg'
  [go]
---------- RUN BEGIN -------------------------------------------------

PARAMETER FILES (1 total):
 0: structure.psf, 570226 atoms, 145437 res, box: Orthorhombic, 144377 mol, 142700 solvent

INPUT TRAJECTORIES (5 total):
 0: 'output.xtc' is a GROMACS XTC file, Parm structure.psf (Orthorhombic box) (reading 100 of 
5000)
 1: 'output.xtc' is a GROMACS XTC file, Parm structure.psf (Orthorhombic box) (reading 100 of 
5000)
 2: 'output.xtc' is a GROMACS XTC file, Parm structure.psf (Orthorhombic box) (reading 100 of 
5000)
 3: 'output.xtc' is a GROMACS XTC file, Parm structure.psf (Orthorhombic box) (reading 100 of 
5000)
 4: 'output.xtc' is a GROMACS XTC file, Parm structure.psf (Orthorhombic box) (reading 100 of 
5000)
  Coordinate processing will occur on 500 frames.

REFERENCE FRAMES (1 total):
    0: MyAvg
        Active reference frame for distance-based masks is 'MyAvg'

PARALLEL INFO:
  Process 0 will handle 32 frames.
  Process 1 will handle 32 frames.
  Process 2 will handle 32 frames.
  Process 3 will handle 32 frames.
  Process 4 will handle 31 frames.
  Process 5 will handle 31 frames.
  Process 6 will handle 31 frames.
  Process 7 will handle 31 frames.
  Process 8 will handle 31 frames.
  Process 9 will handle 31 frames.
  Process 10 will handle 31 frames.
  Process 11 will handle 31 frames.
  Process 12 will handle 31 frames.
  Process 13 will handle 31 frames.
  Process 14 will handle 31 frames.
  Process 15 will handle 31 frames.
.....................................................
ACTION SETUP FOR PARM 'structure.psf' (16 actions):
  0: [strip :SOD,CLA,TIP3,POPC]
        Stripping 553065 atoms.
        Stripped topology: 17161 atoms, 1063 res, box: Orthorhombic, 3 mol
  1: [autoimage]
        Using first molecule as anchor.
        2 molecules are fixed to anchor: 2 3
        0 molecules are mobile.
  2: [rms first @CA]
        Target mask: [@CA](1063)
        Reference topology: structure.psf
        Reference mask: [@CA](1063)
  3: [strip :SOD,CLA,TIP3,POPC]
        Stripping 0 atoms.
Warning: No atoms to strip. Skipping.
Warning: Setup incomplete for [strip]: Skipping
  4: [autoimage]
        Using first molecule as anchor.
        2 molecules are fixed to anchor: 2 3
        0 molecules are mobile.
  5: [rms first @CA]
        Target mask: [@CA](1063)
        Reference topology: structure.psf
        Reference mask: [@CA](1063)
  6: [strip :SOD,CLA,TIP3,POPC]
        Stripping 0 atoms.
Warning: No atoms to strip. Skipping.
Warning: Setup incomplete for [strip]: Skipping
  7: [autoimage]
        Using first molecule as anchor.
        2 molecules are fixed to anchor: 2 3
        0 molecules are mobile.
  8: [rms first @CA]
        Target mask: [@CA](1063)
        Reference topology: structure.psf
        Reference mask: [@CA](1063)
  9: [strip :SOD,CLA,TIP3,POPC]
        Stripping 0 atoms.
Warning: No atoms to strip. Skipping.
Warning: Setup incomplete for [strip]: Skipping
  10: [autoimage]
        Using first molecule as anchor.
        2 molecules are fixed to anchor: 2 3
        0 molecules are mobile.
  11: [rms first @CA]
        Target mask: [@CA](1063)
        Reference topology: structure.psf
        Reference mask: [@CA](1063)
  12: [strip :SOD,CLA,TIP3,POPC]
        Stripping 0 atoms.
Warning: No atoms to strip. Skipping.
Warning: Setup incomplete for [strip]: Skipping
  13: [autoimage]
        Using first molecule as anchor.
        2 molecules are fixed to anchor: 2 3
        0 molecules are mobile.
  14: [rms first @CA]
        Target mask: [@CA](1063)
        Reference topology: structure.psf
        Reference mask: [@CA](1063)
  15: [average crdset MyAvg]
        Mask [*] corresponds to 17161 atoms.
        Averaging over 17161 atoms.

BEGIN PARALLEL TRAJECTORY PROCESSING:
 0% 13% 23% 32% 42% 52% 61% 71% 81% 90% 100% Complete.
TIME: Rank 0 throughput= 25.4529 frames / second.
TIME: Rank 1 throughput= 25.3270 frames / second.
TIME: Rank 2 throughput= 25.1636 frames / second.
TIME: Rank 3 throughput= 25.2121 frames / second.
TIME: Rank 4 throughput= 25.7151 frames / second.
TIME: Rank 5 throughput= 25.5150 frames / second.
TIME: Rank 6 throughput= 25.8017 frames / second.
TIME: Rank 7 throughput= 25.3638 frames / second.
TIME: Rank 8 throughput= 25.5893 frames / second.
TIME: Rank 9 throughput= 25.5936 frames / second.
TIME: Rank 10 throughput= 25.3485 frames / second.
TIME: Rank 11 throughput= 25.1755 frames / second.
TIME: Rank 12 throughput= 25.7129 frames / second.
TIME: Rank 13 throughput= 25.4586 frames / second.
TIME: Rank 14 throughput= 25.5020 frames / second.
TIME: Rank 15 throughput= 25.2399 frames / second.
TIME: Avg. throughput= 393.2017 frames / second.

ACTION OUTPUT:
    AVERAGE: 500 frames, COORDS set 'MyAvg'
Warning: Analysis does not currently use multiple MPI processes.
TIME: Analyses took 0.0001 seconds.

DATASETS (6 total):
        i "i" (string variable), size is 1 (0.033 kB)
        RMSD_00002 "RMSD_00002" (double, rms), size is 500 (4.000 kB)
        RMSD_00003 "RMSD_00003" (double, rms), size is 500 (4.000 kB)
        RMSD_00004 "RMSD_00004" (double, rms), size is 500 (4.000 kB)
        RMSD_00005 "RMSD_00005" (double, rms), size is 500 (4.000 kB)
        RMSD_00006 "RMSD_00006" (double, rms), size is 500 (4.000 kB)
    Total data set memory usage is at least 20.033 kB

RUN TIMING:
TIME:           Init               : 0.2092 s ( 13.36%)
TIME:           Trajectory Process : 1.2716 s ( 81.19%)
TIME:           Data Set Sync      : 0.0012 s (  0.07%)
TIME:           Action Post        : 0.0030 s (  0.19%)
TIME:           Analysis           : 0.0001 s (  0.00%)
TIME:           Data File Write    : 0.0001 s (  0.01%)
TIME:           Other              : 0.0810 s (  0.05%)
TIME:   Run Total 1.5661 s
---------- RUN END ---------------------------------------------------
  [rms ref MyAvg @CA]
Error: Reference mode # atoms on rank (0) != # atoms on master (17161)
Error: Could not initialize action [rms]
        1 errors encountered reading input.
TIME: Total execution time: 6.1552 seconds.
Error: Error(s) occurred during execution.

--------------------------------------------------------------------------
Primary job  terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
--------------------------------------------------------------------------
mpirun detected that one or more processes exited with non-zero status, thus causing
the job to be terminated. The first process to do so was:

  Process name: [[1650,1],10]
  Exit code:    1
--------------------------------------------------------------------------

Why this is not possilbe with mpi? Or i am doing something wrong?

Thanks in advance!

drroe commented 1 month ago

Error: Reference mode # atoms on rank (0) != # atoms on master (17161)

It seems like the calculated average structure isn't getting distributed to all ranks for some reason. I'll open up a pull request to look into this. In the meantime you can split this into two separate runs like so. The first run will be calculating your average structure:

parm prod1/structure.psf
for i=1;i<6;i++
  trajin prod$i/output.xtc 1 5000 50
done
strip :SOD,CLA,TIP3,POPC
autoimage 
rms first @CA
average crdset MyAvg @CA
go
crdout MyAvg myavg.mol2

Note that the strip, autoimage, and rms commands can be outside the loop; they will apply to all incoming trajectories. The loop simply makes it so there are 6 trajin commands, each with a different file name. Once your average is generated you can have a second run to load in the average, RMS-fit to it, and calculate your RMSF values:

parm prod1/structure.psf
parm myavg.mol2
reference myavg.mol2 parm myavg.mol2
for i=1;i<6;i++
  trajin prod$i/output.xtc 1 5000 50
done
strip :SOD,CLA,TIP3,POPC
autoimage 
rms reference @CA
atomicfluct out rmsfA.dat ^3@CA byres
atomicfluct out rmsf465A.dat ^1,@CA byres
run

Both can be run in parallel. Let me know if this makes sense.

H-EKE commented 1 month ago

Hi,

Thanks for your help!! I just tried and it seems to work in 2 different runs.

Also, I had a similar problem with parallel trajout, while using this script

parm prod1/structure.psf
for i=1;i<6;i++
trajin prod$i/output.xtc 1 5000 
strip :SOD,CLA,TIP3,POPC,TIP3,IP3,TIP
autoimage ^1:1-10,^3@CA
rms first @CA
trajout prod$i/cent.xtc
run
clear trajin
done
CPPTRAJ: Trajectory Analysis. V6.4.4 (GitHub) MPI OpenMP
    ___  ___  ___  ___
     | \/ | \/ | \/ | 
    _|_/\_|_/\_|_/\_|_

| Running on 16 processes.
| 1 OpenMP threads available.
| Date/time: 07/22/24 19:29:03
| Available memory: 262.544 GB

INPUT: Reading input from 'STDIN'
  [parm prod1/structure.psf]
        Reading 'prod1/structure.psf' as Charmm PSF
    Reading Charmm PSF file structure.psf as topology file.
        PSF contains 570226 atoms, 145437 residues.
Warning: Determining bond length parameters from element types for 'structure.psf'.
  [for i=1;i<6;i++]
    Setting up 'for' loop.
CONTROL: Parsing control block.
  BLOCK  0: for ((i=1; i<6; i+=1)) do
  [trajin prod$i/output.xtc 1 5000]
        Added command 'trajin' to control block 0.
  [strip :SOD,CLA,TIP3,POPC,TIP3,IP3,TIP]
        Added command 'strip' to control block 0.
  [autoimage ^1:1-10,^3@CA]
        Added command 'autoimage' to control block 0.
  [rms first @CA]
        Added command 'rms' to control block 0.
  [trajout prod$i/cent.xtc]
        Added command 'trajout' to control block 0.
  [run]
        Added command 'run' to control block 0.
  [clear trajin]
        Added command 'clear' to control block 0.
  [done]
  BLOCK  0: END
CONTROL: Executing 1 control block(s).
        Loop over 'i' will execute for 5 iterations.
  [trajin prod1/output.xtc 1 5000]
        Reading 'prod1/output.xtc' as Gromacs XTC
Warning: Trajectory box type is 'Orthorhombic' but topology box type is 'None'.
Warning: Setting topology box information from trajectory.
  [strip :SOD,CLA,TIP3,POPC,TIP3,IP3,TIP]
    STRIP: Stripping atoms in mask [:SOD,CLA,TIP3,POPC,TIP3,IP3,TIP]
  [autoimage ^1:1-10,^3@CA]
    AUTOIMAGE: To box center based on center of mass, anchor mask is [^1:1-10,^3@CA]
  [rms first @CA]
    RMSD: (@CA), reference is first frame (@CA).
        Best-fit RMSD will be calculated, coords will be rotated and translated.
  [trajout prod1/cent.xtc]
        Writing 'prod1/cent.xtc' as Gromacs XTC
  [run]
---------- RUN BEGIN -------------------------------------------------

PARAMETER FILES (1 total):
 0: structure.psf, 570226 atoms, 145437 res, box: Orthorhombic, 144377 mol, 142700 solvent

INPUT TRAJECTORIES (1 total):
 0: 'output.xtc' is a GROMACS XTC file, Parm structure.psf (Orthorhombic box) (reading 5000 of 5257)
  Coordinate processing will occur on 5000 frames.

OUTPUT TRAJECTORIES (1 total):
  'cent.xtc' (5000 frames) is a GROMACS XTC file

PARALLEL INFO:
  Process 0 will handle 313 frames.
  Process 1 will handle 313 frames.
  Process 2 will handle 313 frames.
  Process 3 will handle 313 frames.
  Process 4 will handle 313 frames.
  Process 5 will handle 313 frames.
  Process 6 will handle 313 frames.
  Process 7 will handle 313 frames.
  Process 8 will handle 312 frames.
  Process 9 will handle 312 frames.
  Process 10 will handle 312 frames.
  Process 11 will handle 312 frames.
  Process 12 will handle 312 frames.
  Process 13 will handle 312 frames.
  Process 14 will handle 312 frames.
  Process 15 will handle 312 frames.
.....................................................
ACTION SETUP FOR PARM 'structure.psf' (3 actions):
  0: [strip :SOD,CLA,TIP3,POPC,TIP3,IP3,TIP]
        Stripping 553065 atoms.
        Stripped topology: 17161 atoms, 1063 res, box: Orthorhombic, 3 mol
  1: [autoimage ^1:1-10,^3@CA]
        Anchoring on atoms selected by mask '^1:1-10,^3@CA'
        Mask [^1:1-10,^3@CA] corresponds to 510 atoms.
        3 molecules are fixed to anchor: 1 2 3
        0 molecules are mobile.
  2: [rms first @CA]
        Target mask: [@CA](1063)
        Reference topology: structure.psf
        Reference mask: [@CA](1063)
Error: Could not set up parallel trajout.
Error: Problem with set up or format may be unsupported in parallel.
Error: Setting up output trajectory 'prod1/cent.xtc' in parallel.

RUN TIMING:
TIME:           Init               : 0.0000 s (  0.00%)
TIME:           Trajectory Process : 0.0000 s (  0.00%)
TIME:           Data Set Sync      : 0.0000 s (  0.00%)
TIME:           Action Post        : 0.0000 s (  0.00%)
TIME:           Analysis           : 0.0000 s (  0.00%)
TIME:           Data File Write    : 0.0000 s (  0.00%)
TIME:           Other              : 0.1483 s (  1.00%)
TIME:   Run Total 0.1483 s
---------- RUN END ---------------------------------------------------
        1 errors encountered reading input.
Error: Error(s) occurred during execution.
TIME: Total execution time: 2.5171 seconds.

--------------------------------------------------------------------------
Primary job  terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
--------------------------------------------------------------------------
mpirun detected that one or more processes exited with non-zero status, thus causing
the job to be terminated. The first process to do so was:

  Process name: [[37908,1],0]
  Exit code:    1
--------------------------------------------------------------------------

Is this expected or is it a bug too?

drroe commented 1 month ago

This is expected. Note the earlier error message:

Error: Could not set up parallel trajout.
Error: Problem with set up or format may be unsupported in parallel.
Error: Setting up output trajectory 'prod1/cent.xtc' in parallel.

Trajectory output using the XTC format is not supported in parallel. If you want to output a format in parallel you should use DCD, TRR, NetCDF, etc. If you need a compressed format you can try the NetCDF 4 format with lossy compression (via trajout keywords hdf5 icompress), but this does require compilation with NetCDF4/HDF5. See the manual for full details.

H-EKE commented 1 month ago

Thanks!! :)