Previous Thread
Next Thread
Print Thread
Page 1 of 2 1 2
Image Patch Issue Charmm 36 vs Charmm 34
#31791 04/22/13 06:36 AM
Joined: Sep 2010
Posts: 6
Tigran Offline OP
Forum Member
OP Offline
Forum Member
Joined: Sep 2010
Posts: 6
Dear All,

As I switch from Charmm 34b2 to Charmm 36b1, I notice an error mentioned below. Presumably this has to do with the Image Patch command of Charmm 36b1. Charmm 36b1 is trying to find angle/bond/torsion parameters for the atoms of the polyethylene chain and water molecules. My system consists of 3 layers: top layer - TIP3 water with segid BWAT (6984 atoms), peptide TGTGKGTGT with segid PEP (109 atoms), and a CLA ion with segid CION (1 atom); middle layer - 45 polyethylene chains (4860 atoms); bottom layer - TIP3 water with segid FWAT (4860). Please find CRD, PSF and charmm input files attached to this post in text.tar.gz. Please advise.

Here is the error:

CHARMM> noewald -
CHARMM> ntrfrq 500000
Parameter: DTIMSTEP -> "0.001"
Parameter: DNUMSTEP -> "500000"
IUNREA = -1 IUNWRI = 32 IUNOS = -1
IUNCRD = 32 IUNVEL = -1 KUNIT = -1

SELECTED IMAGES ATOMS BEING CENTERED ABOUT 0.000000 0.000000 0.000000

NONBOND OPTION FLAGS:
ELEC VDW ATOMs CDIElec FSWItch VATOm VFSWIt
BYGRoup NOEXtnd NOEWald
CUTNB = 14.000 CTEXNB =999.000 CTONNB = 8.000 CTOFNB = 12.000
CGONNB = 0.000 CGOFNB = 10.000
WMIN = 1.500 WRNMXD = 0.500 E14FAC = 1.000 EPS = 1.000
NBXMOD = 5
There are 0 atom pairs and 0 atom exclusions.
There are 0 group pairs and 0 group exclusions.
with mode 5 found 24586 exclusions and 14177 interactions(1-4)
found 4662 group exclusions.
: No bond parameters for 1 ( CG321 OT )
: No bond parameters for 2 ( CG321 OT )
: No bond parameters for 3 ( CG321 OT )
: No angle parameters for 1 ( OT CG321 CG321 )
: No angle parameters for 2 ( OT CG321 HGA2 )
: No angle parameters for 3 ( OT CG321 HGA2 )
: No angle parameters for 4 ( OT CG321 CG321 )
: No angle parameters for 5 ( OT CG321 HGA2 )
: No angle parameters for 6 ( OT CG321 HGA2 )
: No angle parameters for 7 ( OT CG321 CG321 )
: No angle parameters for 8 ( OT CG321 HGA2 )
: No angle parameters for 9 ( OT CG321 HGA2 )
: No angle parameters for 136 ( HT OT CG321 )
: No angle parameters for 137 ( HT OT CG321 )
: No angle parameters for 138 ( OT OT CG321 )
: No angle parameters for 139 ( HT OT CG321 )
: No angle parameters for 140 ( HT OT CG321 )
: No angle parameters for 141 ( OT OT CG321 )
: No angle parameters for 142 ( HT OT CG321 )
: No angle parameters for 143 ( HT OT CG321 )
: No angle parameters for 144 ( CG321 OT CG321 )
: No torsion parameters for 1 ( OT CG321 CG321 CG321 )
: No torsion parameters for 2 ( OT OT CG321 CG321 )
: No torsion parameters for 3 ( OT OT OT CG321 )
: No torsion parameters for 4 ( OT CG321 CG321 HGA2 )
: No torsion parameters for 5 ( OT CG321 CG321 HGA2 )
: No torsion parameters for 6 ( OT OT CG321 HGA2 )
: No torsion parameters for 7 ( OT OT CG321 HGA2 )
: No torsion parameters for 8 ( HT OT CG321 CG321 )
: No torsion parameters for 9 ( HT OT CG321 CG321 )
: No torsion parameters for 10 ( HT OT CG321 HGA2 )
: No torsion parameters for 11 ( HT OT CG321 HGA2 )
: No torsion parameters for 12 ( HT OT CG321 HGA2 )
: No torsion parameters for 13 ( HT OT CG321 HGA2 )
: No torsion parameters for 14 ( HT OT OT CG321 )
: No torsion parameters for 15 ( HT OT OT CG321 )
: No torsion parameters for 16 ( OT CG321 CG321 CG321 )
: No torsion parameters for 17 ( OT OT CG321 CG321 )
: No torsion parameters for 18 ( OT OT OT CG321 )
: No torsion parameters for 19 ( OT CG321 CG321 HGA2 )
: No torsion parameters for 20 ( OT CG321 CG321 HGA2 )
: No torsion parameters for 21 ( OT OT CG321 HGA2 )
: No torsion parameters for 22 ( OT OT CG321 HGA2 )
: No torsion parameters for 23 ( HT OT CG321 CG321 )
: No torsion parameters for 24 ( HT OT CG321 CG321 )
: No torsion parameters for 25 ( HT OT CG321 HGA2 )
: And more missing torsion parameters.
: A TOTAL OF 66 MISSING PARAMETERS

***** LEVEL -1 WARNING FROM *****
***** CODES> MISSING PARAMETERS
******************************************
BOMLEV ( 0) IS REACHED - TERMINATING. WRNLEV IS 5



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

And here is my charmm input file:


SET wdir /bioengr/tabramy/polyet/approach1/c34vsc36
SET idir /bioengr/tabramy/polyet/approach1/c34vsc36
SET stdir /bioengr/tabramy/polyet/approach1/c34vsc36
SET odir /bioengr/tabramy/polyet/approach1/c34vsc36

SET rtf1 c22p.rtf
SET rtf2 sams.rtf
SET prm1 c22p.prm
SET prm2 sams.prm

SET psf lys.psf
SET crd lys.crd

SET dtimstep 0.001
SET dnumstep 500000

SET newres dyna1.res
SET newdcd dyna1.dcd
SET newpdb dyna1.pdb
SET newcrd dyna1.crd
SET newpsf dyna1.psf


!********************************************
! Start CHARMM
read rtf card name @wdir/@rtf1
read rtf card append name @wdir/@rtf2
read para card name @wdir/@prm1
read para card append name @wdir/@prm2

read psf card name @idir/@psf
read coor card name @idir/@crd

!********************************************
! Crystallographic indices of the surface
!generation of n cells along U and m cells along V
!specify the overall surface size
set n 18
set m 9
set k 2.5
!unit cell translation parameters
set u1 2.5534
set u2 0.0
set u3 0.0
set v1 0.000000
set v2 4.9491
set v3 0.0
set w1 0.0
set w2 0.0
set w3 7.4241
!Define angle
set theta 90.0000
!Define translational vectors for images
calc bigu1 = @n * @u1
calc bigu2 = @n * @u2
calc bigu3 = @n * @u3
calc bigv1 = @m * @v1
calc bigv2 = @m * @v2
calc bigv3 = @m * @v3
calc bigw1 = @k * @w1
calc bigw2 = @k * @w2
calc bigw3 = @k * @w3


!********************************************
! Box size
SET a0 @bigu1
SET b0 @bigv2
SET c0 68.935

SET theta 90.0
crystal define orthorhombic @a0 @b0 @c0 @theta @theta @theta
crystal build noperations 0 cutoff 15.0
crystal write card unit 6
*

print image trans ! To get the image names


!********************************************
update -
imgfrq -1 -
ixtfrq 0 -
cutim 14.0 -
- ! imal -
- ! NBond-Specifications
inbfrq -1 -
ihbfrq 0 -
- ! bygroup -
update -
- ! Electrostatic Energy Specifications
Elec -
atom -
cdielec -
fswitch -
eps 1.0 -
- ! Van Der Waals Energy Specifications:
VDW -
vatOm -
vfswitch -
- ! Cutoff Distance Specifications:
nbxmod 5 -
e14fac 1.0 -
wmin 1.5 -
ctonnb 8.0 -
ctofnb 12.0 -
cutnb 14.0


!********************************************
! PBC for Polyethylene
update inbf 0 ihbf 0 imgfrq 0 cutim 14.0
! apply patch to chain ends
calc nmax = int (@m * 2 * @k)
set px 0
label loop1
calc px = @px + 1
impa P1A C008 p@px 1 prim p@px 1 setup warn
if @px lt @nmax goto loop1


!********************************************
! Handle Imaging of Atoms
image byres xcen 0.0 ycen 0.0 zcen 0.0 select (resname TIP3 .or. resname PET .or. -
resname CLA) end
image byseg xcen 0.0 ycen 0.0 zcen 0.0 select segid PEP end
image fixed select segid FWAT end


!********************************************
! Setup Constraints for Dynamics Run
cons fix sele segid FWAT end

cons harm mass force 100.0 sele (segid P1 .or. segid P3 .or. segid P5 -
.or. segid P7 .or. segid P9 .or. segid P11 .or. segid P13 .or. segid P15 -
.or. segid P17) .and. .not. type H* end

cons harm mass force 0.5 zscale 0.0 sele (segid P2 .or. segid P4 .or. segid P6 .or. -
segid P8 .or. segid P10 .or. segid P12 .or. segid P14 .or. segid P16 .or. -
segid P18 .or. segid P19 .or. segid P20 .or. segid P21 .or. segid P22 .or. -
segid P23 .or. segid P24 .or. segid P25 .or. segid P26 .or. segid P27 .or. -
segid P28 .or. segid P29 .or. segid P30 .or. segid P31 .or. segid P32 .or. -
segid P33 .or. segid P34 .or. segid P35 .or. segid P36 .or. segid P37 .or. -
segid P38 .or. segid P39 .or. segid P40 .or. segid P41 .or. segid P42 .or. -
segid P43 .or. segid P44 .or. segid P45) .and. type C18 end


!********************************************
! Invoke Shake for MD Run
shake bonh para tol 0.1e-09

open unit 32 write card name @odir/@newres
open unit 32 write file name @odir/@newdcd


!********************************************
! Start Molecular Dynamics Run
dynamics -
vv2 -
timestp @dtimstep -
nstep @dnumstep -
echeck 20.0 -
! firstt 100.0 -
finalt 298.0 -
! teminc 5.0 -
! ihtfrq 2500 -
start -
-! UNIT-SPECifications
! iunrea 30 -
iunwri 32 -
iuncrd 32 -
-! FREQUENCY-SPECifications
inbfrq -1 -
ihbfrq 0 -
imgfrq -1 -
ixtfrq 0 -
iprfrq 1000 -
nprint 1000 -
! nsavc 1000 -
nsavv 1000 -
isvfrq 1000 -
noewald -
ntrfrq 500000


!*******************************************
! Write System Output Files

!close unit 30
close unit 32
close unit 32

coor stat
energy

write coordinates PDB card name @odir/@newpdb
* PDB
*

write coordinates card name @odir/@newcrd
* CRD
*

stop

Attached Files
test.tar.gz (721.85 KB, 186 downloads)
Re: Image Patch Issue Charmm 36 vs Charmm 34
Tigran #31794 04/22/13 07:48 AM
Joined: Sep 2003
Posts: 4,794
Likes: 2
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,794
Likes: 2
Does it work without the impatch command? You have to use the same parameter and topology files as were used when the PSF was generated.

It is easier to follow a shorter input file. Here all you need is the RTF,parameters,PSF, coordinates, crystal setup and a simple "ENERgy" command. Also, avoid using CHARMM variables.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Re: Image Patch Issue Charmm 36 vs Charmm 34
Tigran #31795 04/22/13 03:27 PM
Joined: Oct 2008
Posts: 22
Forum Member
Offline
Forum Member
Joined: Oct 2008
Posts: 22
I work in the same group as the original poster and have been working on the same problem independently.

We had been successfully using the IMPATCH command in c34b2 to create periodic surfaces such as polyethylene and silica. Now that we have switched to c36b1, the same scripts that worked in c34b2 no longer work in c36b1.

For example, I have attached the output for the same polyethylene script run in c34b2 and c36b1. The c34b2 version runs to completion without any noticeable problems while the c36b1 version begins to show problems after the first image update (during minimization) and crashes when it hits dynamics. At this first image update the c36b1 energy spikes from -30,000 kcal/mol (what it should be) to +5 million kcal/mol in the output. When c36b1 reaches dynamics, it complains that it doesn't have the parameters necessary to form covalent bonds between TIP3P waters and the polyethylene surface atoms, which implies something has gone very wrong in c36b1 during that first image update.

The script runs perfectly fine without the IMPATCH command being present, so we have narrowed the problem down to the IMPATCH command in c36b1 failing during the first image update. That, or something is wrong with our input script.

----

Here is the input script used (I can provide other files as necessary):

read rtf card name top_all36_prot_pet.@{parmset}.rtf
read para card name par_all36_prot_pet.@{parmset}.prm
read psf card name petlys.@{parmset}.psf
read coor card name petlys.crd

crystal define orthorhombic 45.9612 44.5419 68.935 90.0 90.0 90.0
crystal build noperations 0 cutoff 15.0
print image trans
crystal write card unit 6
*

image byresidue sele (segid bwat .or. resname pet .or. segid cion) end
image bysegment sele (segid pep) end
image fixed sele (segid fwat) end

update imgfrq -1 ixtfrq 0 cutim 14.0 imall - ! IMAGE SPECs
inbfrq -1 ihbfrq 0 bygroup - ! NBOND SPECs
elec atom cdie fshift eps 1.0 - ! ELEC SPECs
vdw vatom vfswitch - ! VDW SPECs
nbxmod 5 e14fac 1.0 wmin 1.5 - ! CUTOFF SPECs
ctonnb 8.0 ctofnb 12.0 cutnb 14.0

set petchain 1
label loop1
impatch P1A C008 p@petchain 1 prim p@petchain 1 setup warn
incr petchain by 1
if petchain .le. 45 goto loop1

cons fix sele segid FWAT end
cons harm mass force 100.0 sele (segid P1 .or. segid P3 .or. segid P5 -
.or. segid P7 .or. segid P9 .or. segid P11 .or. segid P13 -
.or. segid P15 .or. segid P17) .and. .not. type H* end
cons harm mass force 0.5 zscale 0.0 sele (segid P2 .or. segid P4 .or. segid P6 -
.or. segid P8 .or. segid P10 .or. segid P12 .or. segid P14 .or. segid P16 -
.or. segid P18 .or. segid P19 .or. segid P20 .or. segid P21 .or. segid P22 -
.or. segid P23 .or. segid P24 .or. segid P25 .or. segid P26 .or. segid P27 -
.or. segid P28 .or. segid P29 .or. segid P30 .or. segid P31 .or. segid P32 -
.or. segid P33 .or. segid P34 .or. segid P35 .or. segid P36 .or. segid P37 -
.or. segid P38 .or. segid P39 .or. segid P40 .or. segid P41 .or. segid P42 -
.or. segid P43 .or. segid P44 .or. segid P45) .and. type C18 end

coor stat
energy

mini conj nstep 200
shake bonh param tol 0.1e-09

coor stat
energy

tpcontrol nthermostats 1 ther 1 tref 298.0 tau 0.1 select all end

open unit 31 writ form name @outdir/start.res
open unit 32 writ unform name @outdir/start.dcd

dyna vv2 start iseed 314159 314159 314159 314159 -
nstep 10 timestep 0.001 -
iunread -1 iunwrite 31 iuncrd 32 ntrfrq 9000000 echeck 20 -
iprfrq 10 nprint 10 nsavc 10 nsavv 10 isvfrq 10 -
firstt 48.0 finalt 298.0 teminc 5.0 ihtfrq 100 -
ichecw 1 twindl -5.0 twindh 5.0 iasors 1 iasvel 1 iscvel 0

Attached Files
petoutput.tar.gz (23.24 KB, 182 downloads)
The c36b1 and c34b2 output for the same input script using the IMPATCH command to make a periodic surface. c34b2 runs to completion with no errors while c36b1 begins to show problems after the first image update (during Minimization) and crashes when dyn
Re: Image Patch Issue Charmm 36 vs Charmm 34
Tigran #31796 04/22/13 03:33 PM
Joined: Oct 2008
Posts: 22
Forum Member
Offline
Forum Member
Joined: Oct 2008
Posts: 22
Both Tigran and I used PSFs that are compatible with the topology and parameter files used by the script.

Re: Image Patch Issue Charmm 36 vs Charmm 34
Tigran #31797 04/22/13 06:14 PM
Joined: Sep 2003
Posts: 4,794
Likes: 2
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,794
Likes: 2
Your posts are a bit different. In one case there are clear PSF problems, in the other things go bad after a while in dynamics.
The inputs are differ wrt PSF and IMAGES. All PSF operations (and IMPATCH is a PSF modification) should be done before images are used. As far as I know there have been no changes in the operation of the IMPAtch command. A complex input file is complicated to debug...


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Re: Image Patch Issue Charmm 36 vs Charmm 34
Tigran #31798 04/22/13 06:37 PM
Joined: Sep 2003
Posts: 8,499
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,499
I believe this is an exception, as IMPATCH requires the use of image names, i.e. images must be present prior to the command.


Rick Venable
computational chemist

Re: Image Patch Issue Charmm 36 vs Charmm 34
Tigran #31799 04/23/13 03:31 AM
Joined: Oct 2008
Posts: 22
Forum Member
Offline
Forum Member
Joined: Oct 2008
Posts: 22
Referring to the petoutput.tar.gz that I attached in my first post...

c34b2 runs my script to completion without any crazy nonsense behavior. This gives me some confidence (but not proof) that my script is ok.

c36b1 gives exactly the same output as c34b2 up until the first image update (as vimdiff clearly shows) and then chaos ensues in the c36b1 output. This makes me think that there may be a bug in the c36b1 code in the image updating functions (perhaps an array reading out of bounds???).

Both versions of CHARMM are identical until they diverge in the first image update:

1) c36b1 runs and gives the same number of atoms, groups, and residues in every transformation as c34b2's . So far, so good--identical output.

2) c34b2 then shows:
with mode 5 found 315 exclusions and 675 interactions(1-4)
found 270 image group exclusions.

while c36b1 shows:
with mode 5 found 342 exclusions and 827 interactions(1-4)
found 341 image group exclusions.

The c36b1 energy jumps from a reasonable value of -30,000 kcal/mol to +5 million kcal/mol, going from:

MINI> 100 -29756.21261 158.09396 1.00762 0.95328
MINI INTERN> 1035.04257 829.70962 28.92891 1097.34313 1.50087
MINI CROSS> -21.04765 0.00000 0.00000 0.00000
MINI EXTERN> 843.42907 -31082.83689 0.00000 0.00000 0.00000
MINI IMAGES> -312.95265 -2200.68627 0.00000 0.00000 0.00000
MINI CONSTR> 25.35668 0.00000 0.00000 0.00000 0.00000


to:

MINI> 100 5039908.3146-5069664.5272 2342.8490 0.0000
MINI INTERN> 3270699.4941 1530005.7830 28.9289 1175.3894 1.5009
MINI CROSS> -21.04765 0.00000 0.00000 0.00000
MINI EXTERN> 843.42899 -31082.72123 0.00000 0.00000 0.00000
MINI IMAGES> 270263.29887 -2031.09727 0.00000 0.00000 0.00000
MINI CONSTR> 25.35668 0.00000 0.00000 0.00000 0.00000


during the image update. c34b2 does not exhibit this behavior at all as its energy goes from -29756.21261 (same as c36b1) to -29755.85848 after the update.

If I print out a PSF for my system before and after this image update (not shown in the output I gave), the PSF doesn't change. So, the PSF bonding information isn't altered during the image update. However, there is a large energy jump during the update and when dynamics begins, CHARMM is under the impression that TIP3P waters are suddenly bonded to the polyethylene chain (which is, of course, nonsense) and crashes because the parameter file doesn't have parameters for such covalent bonds.

What can I do to help you determine if this is a bug in c36b1 and not a bug in my complicated script? To me, it's possible that c36b1's MKIMAT2 is not functioning properly with IMPATCH while c34b2's MKIMAT works just fine. The problem could be with MKIMNB or another subroutine as well.

Do I need to show evidence of c36b1 messing up with the IMPATCH command in a simpler system--say one polyethylene chain and water?

Re: Image Patch Issue Charmm 36 vs Charmm 34
Tigran #31800 04/23/13 06:18 AM
Joined: Sep 2003
Posts: 4,794
Likes: 2
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,794
Likes: 2
There is an option which you can try in the dynamics or energy command, IMORiginal. It reverts to using MKIMAT instead of MKIMAT2.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Re: Image Patch Issue Charmm 36 vs Charmm 34
Tigran #31803 04/23/13 07:48 AM
Joined: Oct 2008
Posts: 22
Forum Member
Offline
Forum Member
Joined: Oct 2008
Posts: 22
Thank you for the suggestion, lennart!

When I use IMOR in my ENER, MINI, DYNA, and UPDATE commands, the MKIMAT output agrees perfectly in both c34b2 and c36b1. However, the differences between the c34 and c36 output after the MKIMAT call remain exactly the same.

The entire script output in c36 is identical whether MKIMAT or MKIMAT2 is used (aside from printing the Min-Distance and Upper-Bound values in the last column of the MKIMAT output, of course). MKIMNB in c36 still generates additional nonbonded exclusions and interactions(1-4) that aren't found in c34, c36's MAKGRP finds more image group exclusions than c34, and c36 ends up trying to create the same spurious water-polyethylene bonds when dynamics begins.

It looks as though c36's MKIMAT2 isn't the culprit.

Re: Image Patch Issue Charmm 36 vs Charmm 34
Tigran #31804 04/23/13 07:59 AM
Joined: Sep 2003
Posts: 4,794
Likes: 2
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,794
Likes: 2
Well, it indeed looks like a bug. Please submit a report at www.charmm.org/redmine
I'll see if I can reproduce it with a system of my own.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Re: Image Patch Issue Charmm 36 vs Charmm 34
Tigran #31812 04/23/13 09:04 PM
Joined: Oct 2008
Posts: 22
Forum Member
Offline
Forum Member
Joined: Oct 2008
Posts: 22
Thank you for your help, lennart. I applied for an account at www.charmm.org/redmine and am waiting for it to be activated.

Re: Image Patch Issue Charmm 36 vs Charmm 34
Tigran #31813 04/24/13 08:26 AM
Joined: Oct 2008
Posts: 22
Forum Member
Offline
Forum Member
Joined: Oct 2008
Posts: 22
Tigran (the original poster) and I made a much simpler test case that exhibits the same error using one 'infinite' polyethylene chain in a pure water box.

I used this testcase in the submitted report (issue #129).
Thanks again for your help!

Re: Image Patch Issue Charmm 36 vs Charmm 34
Tigran #31820 04/25/13 01:50 AM
Joined: Sep 2003
Posts: 8,499
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,499
A fix for this problem has been posted in the "Bug Reports and Fixes" forum. CHARMM releases c38b1 and later should not need the fix, but releases c36b1-c37b2 will need it for IMPATCH to work properly.


Rick Venable
computational chemist

Re: Image Patch Issue Charmm 36 vs Charmm 34
Tigran #31821 04/25/13 03:49 AM
Joined: Oct 2008
Posts: 22
Forum Member
Offline
Forum Member
Joined: Oct 2008
Posts: 22
Thank you very much for all your help, Rick and Lennart!

Page 1 of 2 1 2

Moderated by  lennart, rmv 

Link Copied to Clipboard
Powered by UBB.threads™ PHP Forum Software 7.7.4
(Release build 20200307)
Responsive Width:

PHP: 5.6.33-0+deb8u1 Page Time: 0.014s Queries: 44 (0.007s) Memory: 1.0424 MB (Peak: 1.2111 MB) Data Comp: Off Server Time: 2020-09-30 23:03:20 UTC
Valid HTML 5 and Valid CSS