Topic Options
#37257 - 01/03/19 07:55 AM pca data projection problem
dpaul Offline
Forum Member

Registered: 09/05/18
Posts: 2
Loc: india
I am trying to project my pca data in c37b1 version. my input is

* Use QUASI method to generate principal components
* Step 3: generate projections onto selected principal component basis.
* # -Oct-3-1995 (Leo Caves)
*

BOMLev -1
FASTer on

! NB: it is important that you correctly specify the number of modes
! that remain in the mode file. you also need to know how many frames
! are in your trajectory file

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! input
SET system1 x-ind
SET system x
SET avg x-avg
SET sele type CA ! selection
SET nmode 5 ! no. of ("largest") modes kept on file
SET nprj 3 ! no. of ("largest") modes to be projected onto
SET nframe 4000
SET mod bcl2ll-500ns-dif-avg.mod ! mode file
SET infile bcl2ll-500ns-ca-ind.dcd !trajectory
SET average bcl2ll-500ns-ca-ind.crd !average of above (reference for QUASI)

!output
SET prj x-ind.prj
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! setup system
stream toppar.str

! Read PSF
open read unit 10 card name x-ind.psf
read psf unit 10 card

! average set must be in MAIN
OPEN UNIT 1 READ FORM NAME x-ind.crd
READ COOR CARD UNIT 1
CLOSE UNIT 1

! now cycle through trajectory
!OPEN UNIT 21 READ FILE NAME @infile
!TRAJ IREAD 21 NREAD 1

!calculate nframe
!traj quer unit 21
!set nframe ?NFILE

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! for scaling projections
SET natom ?NATOM
! needed for scaling later
SET sqrnat 0
LET sqrnat = SQRT @natom

!no mass weighting
SCALar MASS SET 1.0

OPEN UNIT 77 WRITE FORM NAME @prj
WRITE TITLE UNIT 77
*# principal component projections
*# trajectory file= @infile
*# average coordinates= @average
*# number of frames in trajectory file= @nframe
*# mode file= @mod
*# number of modes to project onto= @nprj (N)
*# id prj1,...,prjN
*

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! load the modes: use TRAVEL's data structures for convenience

SET req 0
SET 1 1
LABEL load

LET req = @nmode - @1
INCR req

OPEN UNIT 21 READ FILE NAME @mod
VIBRAN INBFRQ 0 IHBFRQ 0 NMODE @nmode
NOMASS-weighting
READ NORM UNIT 21
PRINT NORM MODE @req VECTORS ! DOTP
FILL COMP MODE @req FACT 1.0 ! get "largest" mode into COMP set
END

OPEN UNIT 1 WRITE FORM NAME tmp.crd
WRITE COOR CARD UNIT 1 COMP
* junk
*

TREK
TRAJ READ UNIT 25 NOORIENT ! <- the NOOR is IMPORTANT !
tmp.crd
DONE
END

INCR 1
IF 1 LE @nprj GOTO load

OPEN UNIT 1 READ FORM NAME tmp.crd
CLOSE UNIT 1 DISPosition DELETE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! place reference set in COMPARISON
COOR COPY COMP
! keep a copy of them in stores 7-9
SCALar XCOMP STORE 7
SCALar YCOMP STORE 8
SCALar ZCOMP STORE 9

! now cycle through trajectory
OPEN UNIT 25 READ FILE NAME @infile
TRAJ IREAD 25 NREAD 1

!calculate nframe
!traj quer unit 21
!set step ?nstep
!set nframe 0
!let nframe @step / 5000
!


SET 1 1
LABEL LOOP

SET result

! get displacements from reference into MAIN coordinate set
! restore reference set into COMP
SCALAR XCOMP RECALL 7
SCALAR YCOMP RECALL 8
SCALAR ZCOMP RECALL 9

! load frame into MAIN
TRAJ READ
COOR DIFF

SET imod 1
LABEL mode

TRAVEL ! put current mode in COMP coordinate set
COPY INDE @imod COMP
END

! load displacements into the stores
SCALAR X STORE 1
SCALAR Y STORE 2
SCALAR Z STORE 3
! get product of displacements and modes
SCALAR XCOMP *STORE 1
SCALAR YCOMP *STORE 2
SCALAR ZCOMP *STORE 3
! sum up the 3 components (FBETA as scratch array!)
SCALAR FBETA RECALL 2
SCALAR FBETA +STORE 1
SCALAR FBETA RECALL 3
SCALAR FBETA +STORE 1
! finally unload result and get sum
SCALAR FBETA RECALL 1
SCALAR FBETA STATistics
! we now have the dot product ! scale such that the distance reflects
! an RMSD of structures. therefore we want to divide by sqrt(natoms)
SET dot ?STOT
DIVI dot BY @sqrnat

SET result @result @dot

INCR imod
IF imod LE @nprj GOTO mode

WRITE TITLE UNIT 77
* @1 @result
*

INCR 1
IF 1 LE @nframe GOTO LOOP

STOP

while running this I am getting the error

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0 0x7F5BEE426697
#1 0x7F5BEE426CDE
#2 0x7F5BED92D24F
#3 0xE8E197 in __travelmain_MOD_trekcom
#4 0xE986E8 in __travelmain_MOD_trek
#5 0x414372 in maincomx_
#6 0x416D54 in MAIN__ at charmm_main.src:?
Segmentation fault (core dumped)

Please help me. Thanks in advance.

Top
#37258 - 01/03/19 11:55 AM Re: pca data projection problem [Re: dpaul]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8336
Loc: 39 03 48 N, 77 06 54 W
Segmentation faults can arise from many causes, and can be merely the final gasp of an earlier problem. Careful reading of the output log file for other warning or error messages should be your next step.

The problem could be bad input of some kind (syntax, numerical value, data file, etc.) or it could be a bug, fixed in a later CHARMM version; c43b1 is available, and the free version may meet your needs . There can be bugs that are compiler dependent, so it's often worth trying a different compiler.

_________________________
Rick Venable
computational chemist


Top

Moderator:  BRBrooks, John Legato