8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors
534 write (iout,*) 'FTORS factor =',ftors
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
556 if ( secstruc(i) .eq. 'H') then
557 ! Helix restraints for this residue
560 phi0(ii) = 45.0D0*deg2rad
561 drange(ii)= 5.0D0*deg2rad
562 phibound(1,i) = phi0(ii)-drange(ii)
563 phibound(2,i) = phi0(ii)+drange(ii)
564 else if (secstruc(i) .eq. 'E') then
565 ! strand restraints for this residue
568 phi0(ii) = 180.0D0*deg2rad
569 drange(ii)= 5.0D0*deg2rad
570 phibound(1,i) = phi0(ii)-drange(ii)
571 phibound(2,i) = phi0(ii)+drange(ii)
573 ! no restraints for this residue
574 ndih_nconstr=ndih_nconstr+1
575 idih_nconstr(ndih_nconstr)=i
579 ! deallocate(secstruc)
582 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
583 ! deallocate(secstruc)
586 write(iout,'(A20)')'Error reading FTORS'
587 ! deallocate(secstruc)
589 end subroutine secstrp2dihc
590 !-----------------------------------------------------------------------------
591 subroutine read_secstr_pred(jin,jout,errors)
593 ! implicit real*8 (a-h,o-z)
594 ! INCLUDE 'DIMENSIONS'
595 ! include 'COMMON.IOUNITS'
596 ! include 'COMMON.CHAIN'
597 !el character(len=1),dimension(nres) :: secstruc !(maxres)
598 !el COMMON/SECONDARYS/secstruc
600 character(len=80) :: line,line1 !,ucase
601 logical :: errflag,errors,blankline
604 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
607 read (jin,'(a)') line
608 write (jout,'(2a)') '> ',line(1:78)
610 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
611 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
612 ! end-groups in the input file "*.spred"
615 do while (index(line1,'$END').eq.0)
616 ! Override commented lines.
619 do while (.not.blankline)
621 call mykey(line,line1,ipos,blankline,errflag)
622 if (errflag) write (jout,'(2a)') &
623 'Error when reading sequence in line: ',line
624 errors=errors .or. errflag
625 if (.not. blankline .and. .not. errflag) then
628 !el if (iseq.le.maxres) then
629 if (line1(1:1).eq.'-' ) then
630 secstruc(iseq)=line1(1:1)
631 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
632 ( ucase(line1(1:1)).eq.'H' ) ) then
633 secstruc(iseq)=ucase(line1(1:1))
636 write (jout,1010) line1(1:1), iseq
641 !el write (jout,1000) iseq,maxres
644 do while (ipos1.le.iend)
649 !el if (iseq.le.maxres) then
650 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
651 secstruc(iseq)=line1(ipos1-1:ipos1-1)
652 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
653 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
654 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
657 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
662 !el write (jout,1000) iseq,maxres
669 read (jin,'(a)') line
670 write (jout,'(2a)') '> ',line(1:78)
674 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
676 !d check whether the found length of the chain is correct.
677 length_of_chain=iseq-1
678 if (length_of_chain .ne. nres) then
680 write (jout,'(a,i4,a,i4,a)') &
681 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
682 ,length_of_chain,') does not match with the number of residues (' &
687 1000 format('Error - the number of residues (',i4,&
688 ') has exceeded maximum (',i4,').')
689 1010 format ('Error - unrecognized secondary structure label',a4,&
692 end subroutine read_secstr_pred
694 !-----------------------------------------------------------------------------
696 !-----------------------------------------------------------------------------
701 use control_data, only:maxterm !,maxtor
705 use control, only: getenv_loc
707 ! Read the parameters of the probability distributions of the virtual-bond
708 ! valence angles and the side chains and energy parameters.
710 ! Important! Energy-term weights ARE NOT read here; they are read from the
711 ! main input file instead, because NO defaults have yet been set for these
714 ! implicit real*8 (a-h,o-z)
715 ! include 'DIMENSIONS'
720 ! include 'COMMON.IOUNITS'
721 ! include 'COMMON.CHAIN'
722 ! include 'COMMON.INTERACT'
723 ! include 'COMMON.GEO'
724 ! include 'COMMON.LOCAL'
725 ! include 'COMMON.TORSION'
726 ! include 'COMMON.SCCOR'
727 ! include 'COMMON.SCROT'
728 ! include 'COMMON.FFIELD'
729 ! include 'COMMON.NAMES'
730 ! include 'COMMON.SBRIDGE'
731 ! include 'COMMON.MD'
732 ! include 'COMMON.SETUP'
733 character(len=1) :: t1,t2,t3
734 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
735 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
736 logical :: lprint,LaTeX
737 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
738 real(kind=8),dimension(13) :: b
739 character(len=3) :: lancuch !,ucase
741 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
742 integer :: maxinter,junk,kk,ii
743 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
744 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
745 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
746 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
748 integer :: ichir1,ichir2
749 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
750 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
751 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
753 ! For printing parameters after they are read set the following in the UNRES
756 ! setenv PRINT_PARM YES
758 ! To print parameters in LaTeX format rather than as ASCII tables:
762 call getenv_loc("PRINT_PARM",lancuch)
763 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
764 call getenv_loc("LATEX",lancuch)
765 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 dwa16=2.0d0**(1.0d0/6.0d0)
769 ! Assign virtual-bond length
772 vblinv2=vblinv*vblinv
774 ! Read the virtual-bond parameters, masses, and moments of inertia
775 ! and Stokes' radii of the peptide group and side chains
777 allocate(dsc(ntyp1)) !(ntyp1)
778 allocate(dsc_inv(ntyp1)) !(ntyp1)
779 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
780 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
781 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
782 allocate(nbondterm(ntyp)) !(ntyp)
783 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
784 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
785 allocate(msc(ntyp+1,5)) !(ntyp+1)
786 allocate(isc(ntyp+1,5)) !(ntyp+1)
787 allocate(restok(ntyp+1,5)) !(ntyp+1)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 read (ibond,*) vbldp0,akp,mp,ip,pstok
798 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
799 dsc(i) = vbldsc0(1,i)
803 dsc_inv(i)=1.0D0/dsc(i)
807 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
809 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
810 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
811 dsc(i) = vbldsc0(1,i)
815 dsc_inv(i)=1.0D0/dsc(i)
819 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
822 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
823 ! dsc(i) = vbldsc0_nucl(1,i)
827 ! dsc_inv(i)=1.0D0/dsc(i)
830 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
831 ! do i=1,ntyp_molec(2)
832 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
833 ! aksc_nucl(j,i),abond0_nucl(j,i),&
834 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
835 ! dsc(i) = vbldsc0(1,i)
839 ! dsc_inv(i)=1.0D0/dsc(i)
844 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
845 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
847 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
849 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
850 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
852 write (iout,'(13x,3f10.5)') &
853 vbldsc0(j,i),aksc(j,i),abond0(j,i)
857 !----------------------------------------------------
858 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
859 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
860 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
861 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
862 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
863 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
873 allocate(liptranene(ntyp))
874 !C reading lipid parameters
875 write (iout,*) "iliptranpar",iliptranpar
877 read(iliptranpar,*) pepliptran
880 read(iliptranpar,*) liptranene(i)
881 print *,liptranene(i)
887 ! Read the parameters of the probability distribution/energy expression
888 ! of the virtual-bond valence angles theta
891 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
892 (bthet(j,i,1,1),j=1,2)
893 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
894 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
895 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
899 athet(1,i,1,-1)=athet(1,i,1,1)
900 athet(2,i,1,-1)=athet(2,i,1,1)
901 bthet(1,i,1,-1)=-bthet(1,i,1,1)
902 bthet(2,i,1,-1)=-bthet(2,i,1,1)
903 athet(1,i,-1,1)=-athet(1,i,1,1)
904 athet(2,i,-1,1)=-athet(2,i,1,1)
905 bthet(1,i,-1,1)=bthet(1,i,1,1)
906 bthet(2,i,-1,1)=bthet(2,i,1,1)
910 athet(1,i,-1,-1)=athet(1,-i,1,1)
911 athet(2,i,-1,-1)=-athet(2,-i,1,1)
912 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
913 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
914 athet(1,i,-1,1)=athet(1,-i,1,1)
915 athet(2,i,-1,1)=-athet(2,-i,1,1)
916 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
917 bthet(2,i,-1,1)=bthet(2,-i,1,1)
918 athet(1,i,1,-1)=-athet(1,-i,1,1)
919 athet(2,i,1,-1)=athet(2,-i,1,1)
920 bthet(1,i,1,-1)=bthet(1,-i,1,1)
921 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
926 polthet(j,i)=polthet(j,-i)
929 gthet(j,i)=gthet(j,-i)
937 'Parameters of the virtual-bond valence angles:'
938 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
939 ' ATHETA0 ',' A1 ',' A2 ',&
942 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
943 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
945 write (iout,'(/a/9x,5a/79(1h-))') &
946 'Parameters of the expression for sigma(theta_c):',&
947 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
948 ' ALPH3 ',' SIGMA0C '
950 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
951 (polthet(j,i),j=0,3),sigc0(i)
953 write (iout,'(/a/9x,5a/79(1h-))') &
954 'Parameters of the second gaussian:',&
955 ' THETA0 ',' SIGMA0 ',' G1 ',&
958 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
959 sig0(i),(gthet(j,i),j=1,3)
963 'Parameters of the virtual-bond valence angles:'
964 write (iout,'(/a/9x,5a/79(1h-))') &
965 'Coefficients of expansion',&
966 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
967 ' b1*10^1 ',' b2*10^1 '
969 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
970 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
971 (10*bthet(j,i,1,1),j=1,2)
973 write (iout,'(/a/9x,5a/79(1h-))') &
974 'Parameters of the expression for sigma(theta_c):',&
975 ' alpha0 ',' alph1 ',' alph2 ',&
976 ' alhp3 ',' sigma0c '
978 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
979 (polthet(j,i),j=0,3),sigc0(i)
981 write (iout,'(/a/9x,5a/79(1h-))') &
982 'Parameters of the second gaussian:',&
983 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
986 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
987 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
993 ! Read the parameters of Utheta determined from ab initio surfaces
994 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
996 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
997 ntheterm3,nsingle,ndouble
998 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1000 !----------------------------------------------------
1001 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1002 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1003 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1004 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1005 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1006 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1007 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1008 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1009 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1010 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1011 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1012 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1013 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1014 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1015 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1016 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1017 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1018 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1019 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1020 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1021 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1022 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1023 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1024 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1026 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1028 ithetyp(i)=-ithetyp(-i)
1031 aa0thet(:,:,:,:)=0.0d0
1032 aathet(:,:,:,:,:)=0.0d0
1033 bbthet(:,:,:,:,:,:)=0.0d0
1034 ccthet(:,:,:,:,:,:)=0.0d0
1035 ddthet(:,:,:,:,:,:)=0.0d0
1036 eethet(:,:,:,:,:,:)=0.0d0
1037 ffthet(:,:,:,:,:,:,:)=0.0d0
1038 ggthet(:,:,:,:,:,:,:)=0.0d0
1040 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1042 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1043 ! VAR:1=non-glicyne non-proline 2=proline
1044 ! VAR:negative values for D-aminoacid
1046 do j=-nthetyp,nthetyp
1047 do k=-nthetyp,nthetyp
1048 read (ithep,'(6a)',end=111,err=111) res1
1049 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1050 ! VAR: aa0thet is variable describing the average value of Foureir
1051 ! VAR: expansion series
1052 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1053 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1054 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1055 read (ithep,*,end=111,err=111) &
1056 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1057 read (ithep,*,end=111,err=111) &
1058 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1059 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1060 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1061 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1063 read (ithep,*,end=111,err=111) &
1064 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1065 ffthet(lll,llll,ll,i,j,k,iblock),&
1066 ggthet(llll,lll,ll,i,j,k,iblock),&
1067 ggthet(lll,llll,ll,i,j,k,iblock),&
1068 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1073 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1074 ! coefficients of theta-and-gamma-dependent terms are zero.
1075 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1076 ! RECOMENTDED AFTER VERSION 3.3)
1080 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1081 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1083 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1084 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1087 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1089 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1092 ! AND COMMENT THE LOOPS BELOW
1096 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1097 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1099 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1100 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1103 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1105 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1110 ! Substitution for D aminoacids from symmetry.
1113 do j=-nthetyp,nthetyp
1114 do k=-nthetyp,nthetyp
1115 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1117 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1121 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1122 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1123 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1124 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1130 ffthet(llll,lll,ll,i,j,k,iblock)= &
1131 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1132 ffthet(lll,llll,ll,i,j,k,iblock)= &
1133 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1134 ggthet(llll,lll,ll,i,j,k,iblock)= &
1135 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1136 ggthet(lll,llll,ll,i,j,k,iblock)= &
1137 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1146 ! Control printout of the coefficients of virtual-bond-angle potentials
1149 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1154 write (iout,'(//4a)') &
1155 'Type ',onelett(i),onelett(j),onelett(k)
1156 write (iout,'(//a,10x,a)') " l","a[l]"
1157 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1158 write (iout,'(i2,1pe15.5)') &
1159 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1161 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1162 "b",l,"c",l,"d",l,"e",l
1164 write (iout,'(i2,4(1pe15.5))') m,&
1165 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1166 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1170 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1171 "f+",l,"f-",l,"g+",l,"g-",l
1174 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1175 ffthet(n,m,l,i,j,k,iblock),&
1176 ffthet(m,n,l,i,j,k,iblock),&
1177 ggthet(n,m,l,i,j,k,iblock),&
1178 ggthet(m,n,l,i,j,k,iblock)
1188 write (2,*) "Start reading THETA_PDB",ithep_pdb
1190 ! write (2,*) 'i=',i
1191 read (ithep_pdb,*,err=111,end=111) &
1192 a0thet(i),(athet(j,i,1,1),j=1,2),&
1193 (bthet(j,i,1,1),j=1,2)
1194 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1195 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1196 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1197 sigc0(i)=sigc0(i)**2
1200 athet(1,i,1,-1)=athet(1,i,1,1)
1201 athet(2,i,1,-1)=athet(2,i,1,1)
1202 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1203 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1204 athet(1,i,-1,1)=-athet(1,i,1,1)
1205 athet(2,i,-1,1)=-athet(2,i,1,1)
1206 bthet(1,i,-1,1)=bthet(1,i,1,1)
1207 bthet(2,i,-1,1)=bthet(2,i,1,1)
1210 a0thet(i)=a0thet(-i)
1211 athet(1,i,-1,-1)=athet(1,-i,1,1)
1212 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1213 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1214 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1215 athet(1,i,-1,1)=athet(1,-i,1,1)
1216 athet(2,i,-1,1)=-athet(2,-i,1,1)
1217 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1218 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1219 athet(1,i,1,-1)=-athet(1,-i,1,1)
1220 athet(2,i,1,-1)=athet(2,-i,1,1)
1221 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1222 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1223 theta0(i)=theta0(-i)
1227 polthet(j,i)=polthet(j,-i)
1230 gthet(j,i)=gthet(j,-i)
1233 write (2,*) "End reading THETA_PDB"
1237 !--------------- Reading theta parameters for nucleic acid-------
1238 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1239 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1240 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1241 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1242 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1243 nthetyp_nucl+1,nthetyp_nucl+1))
1244 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1245 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1246 nthetyp_nucl+1,nthetyp_nucl+1))
1247 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1248 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1249 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1250 nthetyp_nucl+1,nthetyp_nucl+1))
1251 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1252 nthetyp_nucl+1,nthetyp_nucl+1))
1253 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1254 nthetyp_nucl+1,nthetyp_nucl+1))
1255 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1256 nthetyp_nucl+1,nthetyp_nucl+1))
1257 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1258 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1259 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1260 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1261 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1262 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1264 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1265 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1267 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1269 aa0thet_nucl(:,:,:)=0.0d0
1270 aathet_nucl(:,:,:,:)=0.0d0
1271 bbthet_nucl(:,:,:,:,:)=0.0d0
1272 ccthet_nucl(:,:,:,:,:)=0.0d0
1273 ddthet_nucl(:,:,:,:,:)=0.0d0
1274 eethet_nucl(:,:,:,:,:)=0.0d0
1275 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1276 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1281 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1282 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1283 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1284 read (ithep_nucl,*,end=111,err=111) &
1285 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1286 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1287 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1288 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1289 read (ithep_nucl,*,end=111,err=111) &
1290 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1291 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1292 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1297 !-------------------------------------------
1298 allocate(nlob(ntyp1)) !(ntyp1)
1299 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1300 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1301 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1308 gaussc(:,:,:,:)=0.0D0
1312 ! Read the parameters of the probability distribution/energy expression
1313 ! of the side chains.
1316 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1320 dsc_inv(i)=1.0D0/dsc(i)
1331 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1332 ((blower(k,l,1),l=1,k),k=1,3)
1333 censc(1,1,-i)=censc(1,1,i)
1334 censc(2,1,-i)=censc(2,1,i)
1335 censc(3,1,-i)=-censc(3,1,i)
1337 read (irotam,*,end=112,err=112) bsc(j,i)
1338 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1339 ((blower(k,l,j),l=1,k),k=1,3)
1340 censc(1,j,-i)=censc(1,j,i)
1341 censc(2,j,-i)=censc(2,j,i)
1342 censc(3,j,-i)=-censc(3,j,i)
1343 ! BSC is amplitude of Gaussian
1350 akl=akl+blower(k,m,j)*blower(l,m,j)
1354 if (((k.eq.3).and.(l.ne.3)) &
1355 .or.((l.eq.3).and.(k.ne.3))) then
1356 gaussc(k,l,j,-i)=-akl
1357 gaussc(l,k,j,-i)=-akl
1359 gaussc(k,l,j,-i)=akl
1360 gaussc(l,k,j,-i)=akl
1369 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1372 if (nlobi.gt.0) then
1374 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1375 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1376 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1377 'log h',(bsc(j,i),j=1,nlobi)
1378 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1379 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1381 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1382 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1385 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1386 write (iout,'(a,f10.4,4(16x,f10.4))') &
1387 'Center ',(bsc(j,i),j=1,nlobi)
1388 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1397 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1398 ! added by Urszula Kozlowska 07/11/2007
1400 !el Maximum number of SC local term fitting function coefficiants
1401 !el integer,parameter :: maxsccoef=65
1403 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1406 read (irotam,*,end=112,err=112)
1408 read (irotam,*,end=112,err=112)
1411 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1415 !---------reading nucleic acid parameters for rotamers-------------------
1416 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1417 do i=1,ntyp_molec(2)
1418 read (irotam_nucl,*,end=112,err=112)
1420 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1426 write (iout,*) "Base rotamer parameters"
1427 do i=1,ntyp_molec(2)
1428 write (iout,'(a)') restyp(i,2)
1429 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1434 ! Read the parameters of the probability distribution/energy expression
1435 ! of the side chains.
1437 write (2,*) "Start reading ROTAM_PDB"
1439 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1443 dsc_inv(i)=1.0D0/dsc(i)
1454 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1455 ((blower(k,l,1),l=1,k),k=1,3)
1457 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1458 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1459 ((blower(k,l,j),l=1,k),k=1,3)
1466 akl=akl+blower(k,m,j)*blower(l,m,j)
1476 write (2,*) "End reading ROTAM_PDB"
1482 ! Read torsional parameters in old format
1484 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1486 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1487 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1488 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1490 !el from energy module--------
1491 allocate(v1(nterm_old,ntortyp,ntortyp))
1492 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1493 !el---------------------------
1498 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1504 write (iout,'(/a/)') 'Torsional constants:'
1507 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1508 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1514 ! Read torsional parameters
1516 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1517 read (itorp,*,end=113,err=113) ntortyp
1518 !el from energy module---------
1519 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1520 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1522 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1523 allocate(vlor2(maxlor,ntortyp,ntortyp))
1524 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1525 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1527 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1528 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1529 !el---------------------------
1532 !el---------------------------
1534 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1536 itortyp(i)=-itortyp(-i)
1538 itortyp(ntyp1)=ntortyp
1539 itortyp(-ntyp1)=-ntortyp
1541 write (iout,*) 'ntortyp',ntortyp
1543 do j=-ntortyp+1,ntortyp-1
1544 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1546 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1547 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1550 do k=1,nterm(i,j,iblock)
1551 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1553 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1554 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1555 v0ij=v0ij+si*v1(k,i,j,iblock)
1557 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1558 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1559 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1561 do k=1,nlor(i,j,iblock)
1562 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1563 vlor2(k,i,j),vlor3(k,i,j)
1564 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1567 v0(-i,-j,iblock)=v0ij
1573 write (iout,'(/a/)') 'Torsional constants:'
1575 do i=-ntortyp,ntortyp
1576 do j=-ntortyp,ntortyp
1577 write (iout,*) 'ityp',i,' jtyp',j
1578 write (iout,*) 'Fourier constants'
1579 do k=1,nterm(i,j,iblock)
1580 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1583 write (iout,*) 'Lorenz constants'
1584 do k=1,nlor(i,j,iblock)
1585 write (iout,'(3(1pe15.5))') &
1586 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1592 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1594 ! 6/23/01 Read parameters for double torsionals
1596 !el from energy module------------
1597 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1598 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1599 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1600 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1601 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1602 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1603 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1604 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1605 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1606 !---------------------------------
1610 do j=-ntortyp+1,ntortyp-1
1611 do k=-ntortyp+1,ntortyp-1
1612 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1613 ! write (iout,*) "OK onelett",
1616 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1617 .or. t3.ne.toronelet(k)) then
1618 write (iout,*) "Error in double torsional parameter file",&
1621 call MPI_Finalize(Ierror)
1623 stop "Error in double torsional parameter file"
1625 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1626 ntermd_2(i,j,k,iblock)
1627 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1628 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1629 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1630 ntermd_1(i,j,k,iblock))
1631 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1632 ntermd_1(i,j,k,iblock))
1633 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1634 ntermd_1(i,j,k,iblock))
1635 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1636 ntermd_1(i,j,k,iblock))
1637 ! Martix of D parameters for one dimesional foureir series
1638 do l=1,ntermd_1(i,j,k,iblock)
1639 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1640 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1641 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1642 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1643 ! write(iout,*) "whcodze" ,
1644 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1646 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1647 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1648 v2s(m,l,i,j,k,iblock),&
1649 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1650 ! Martix of D parameters for two dimesional fourier series
1651 do l=1,ntermd_2(i,j,k,iblock)
1653 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1654 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1655 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1656 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1665 write (iout,*) 'Constants for double torsionals'
1668 do j=-ntortyp+1,ntortyp-1
1669 do k=-ntortyp+1,ntortyp-1
1670 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1671 ' nsingle',ntermd_1(i,j,k,iblock),&
1672 ' ndouble',ntermd_2(i,j,k,iblock)
1674 write (iout,*) 'Single angles:'
1675 do l=1,ntermd_1(i,j,k,iblock)
1676 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1677 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1678 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1679 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1682 write (iout,*) 'Pairs of angles:'
1683 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1684 do l=1,ntermd_2(i,j,k,iblock)
1685 write (iout,'(i5,20f10.5)') &
1686 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1689 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1690 do l=1,ntermd_2(i,j,k,iblock)
1691 write (iout,'(i5,20f10.5)') &
1692 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1693 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1702 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1703 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
1704 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
1705 !el from energy module---------
1706 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1707 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1709 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1710 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
1711 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
1712 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1714 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
1715 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1716 !el---------------------------
1719 !el--------------------
1720 read (itorp_nucl,*,end=113,err=113) &
1721 (itortyp_nucl(i),i=1,ntyp_molec(2))
1722 ! print *,itortyp_nucl(:)
1723 !c write (iout,*) 'ntortyp',ntortyp
1726 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
1727 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
1730 do k=1,nterm_nucl(i,j)
1731 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
1732 v0ij=v0ij+si*v1_nucl(k,i,j)
1735 do k=1,nlor_nucl(i,j)
1736 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
1737 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
1738 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
1744 ! Read of Side-chain backbone correlation parameters
1745 ! Modified 11 May 2012 by Adasko
1748 read (isccor,*,end=119,err=119) nsccortyp
1750 !el from module energy-------------
1751 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1752 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1753 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1754 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1755 !-----------------------------------
1757 !el from module energy-------------
1758 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1760 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1762 isccortyp(i)=-isccortyp(-i)
1764 iscprol=isccortyp(20)
1765 ! write (iout,*) 'ntortyp',ntortyp
1767 !c maxinter is maximum interaction sites
1768 !el from module energy---------
1769 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1770 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1771 -nsccortyp:nsccortyp))
1772 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1773 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1774 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1775 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1776 !-----------------------------------
1780 read (isccor,*,end=119,err=119) &
1781 nterm_sccor(i,j),nlor_sccor(i,j)
1787 nterm_sccor(-i,j)=nterm_sccor(i,j)
1788 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1789 nterm_sccor(i,-j)=nterm_sccor(i,j)
1790 do k=1,nterm_sccor(i,j)
1791 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1793 if (j.eq.iscprol) then
1794 if (i.eq.isccortyp(10)) then
1795 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1796 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1798 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1799 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1800 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1801 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1802 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1803 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1804 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1805 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1808 if (i.eq.isccortyp(10)) then
1809 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1810 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1812 if (j.eq.isccortyp(10)) then
1813 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1814 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1816 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1817 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1818 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1819 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1820 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1821 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1825 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1826 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1827 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1828 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1831 do k=1,nlor_sccor(i,j)
1832 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1833 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1834 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1835 (1+vlor3sccor(k,i,j)**2)
1837 v0sccor(l,i,j)=v0ijsccor
1838 v0sccor(l,-i,j)=v0ijsccor1
1839 v0sccor(l,i,-j)=v0ijsccor2
1840 v0sccor(l,-i,-j)=v0ijsccor3
1846 !el from module energy-------------
1847 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1849 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1850 ! write (iout,*) 'ntortyp',ntortyp
1852 !c maxinter is maximum interaction sites
1853 !el from module energy---------
1854 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1855 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1856 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1857 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1858 !-----------------------------------
1862 read (isccor,*,end=119,err=119) &
1863 nterm_sccor(i,j),nlor_sccor(i,j)
1867 do k=1,nterm_sccor(i,j)
1868 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1870 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1873 do k=1,nlor_sccor(i,j)
1874 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1875 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1876 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1877 (1+vlor3sccor(k,i,j)**2)
1879 v0sccor(l,i,j)=v0ijsccor !el ,iblock
1887 write (iout,'(/a/)') 'Torsional constants:'
1890 write (iout,*) 'ityp',i,' jtyp',j
1891 write (iout,*) 'Fourier constants'
1892 do k=1,nterm_sccor(i,j)
1893 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1895 write (iout,*) 'Lorenz constants'
1896 do k=1,nlor_sccor(i,j)
1897 write (iout,'(3(1pe15.5))') &
1898 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1905 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1906 ! interaction energy of the Gly, Ala, and Pro prototypes.
1911 write (iout,*) "Coefficients of the cumulants"
1913 read (ifourier,*) nloctyp
1914 !write(iout,*) "nloctyp",nloctyp
1915 !el from module energy-------
1916 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1917 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1918 allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1919 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1920 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1921 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1922 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1923 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1927 !--------------------------------
1930 read (ifourier,*,end=115,err=115)
1931 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
1933 write (iout,*) 'Type',i
1934 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
1944 B1tilde(1,-i) =-b(3)
1946 ! b1tilde(1,i)=0.0d0
1947 ! b1tilde(2,i)=0.0d0
1972 Ctilde(1,2,-i)=-b(9)
1976 ! Ctilde(1,1,i)=0.0d0
1977 ! Ctilde(1,2,i)=0.0d0
1978 ! Ctilde(2,1,i)=0.0d0
1979 ! Ctilde(2,2,i)=0.0d0
1997 Dtilde(1,2,-i)=-b(8)
2001 ! Dtilde(1,1,i)=0.0d0
2002 ! Dtilde(1,2,i)=0.0d0
2003 ! Dtilde(2,1,i)=0.0d0
2004 ! Dtilde(2,2,i)=0.0d0
2005 EE(1,1,i)= b(10)+b(11)
2006 EE(2,2,i)=-b(10)+b(11)
2007 EE(2,1,i)= b(12)-b(13)
2008 EE(1,2,i)= b(12)+b(13)
2009 EE(1,1,-i)= b(10)+b(11)
2010 EE(2,2,-i)=-b(10)+b(11)
2011 EE(2,1,-i)=-b(12)+b(13)
2012 EE(1,2,-i)=-b(12)-b(13)
2018 ! ee(2,1,i)=ee(1,2,i)
2022 write (iout,*) 'Type',i
2024 write(iout,*) B1(1,i),B1(2,i)
2026 write(iout,*) B2(1,i),B2(2,i)
2029 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
2033 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
2037 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
2042 ! Read electrostatic-interaction parameters
2047 write (iout,'(/a)') 'Electrostatic interaction constants:'
2048 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2049 'IT','JT','APP','BPP','AEL6','AEL3'
2051 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2052 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2053 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2054 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2059 app (i,j)=epp(i,j)*rri*rri
2060 bpp (i,j)=-2.0D0*epp(i,j)*rri
2061 ael6(i,j)=elpp6(i,j)*4.2D0**6
2062 ael3(i,j)=elpp3(i,j)*4.2D0**3
2064 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2070 ! Read side-chain interaction parameters.
2072 !el from module energy - COMMON.INTERACT-------
2073 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2074 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2075 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2077 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2078 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2079 allocate(epslip(ntyp,ntyp))
2087 !--------------------------------
2089 read (isidep,*,end=117,err=117) ipot,expon
2090 if (ipot.lt.1 .or. ipot.gt.5) then
2091 write (iout,'(2a)') 'Error while reading SC interaction',&
2092 'potential file - unknown potential type.'
2094 call MPI_Finalize(Ierror)
2100 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2101 ', exponents are ',expon,2*expon
2102 ! goto (10,20,30,30,40) ipot
2104 !----------------------- LJ potential ---------------------------------
2106 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2107 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2108 (sigma0(i),i=1,ntyp)
2110 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2111 write (iout,'(a/)') 'The epsilon array:'
2112 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2113 write (iout,'(/a)') 'One-body parameters:'
2114 write (iout,'(a,4x,a)') 'residue','sigma'
2115 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2118 !----------------------- LJK potential --------------------------------
2120 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2121 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2122 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2124 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2125 write (iout,'(a/)') 'The epsilon array:'
2126 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2127 write (iout,'(/a)') 'One-body parameters:'
2128 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2129 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2133 !---------------------- GB or BP potential -----------------------------
2137 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2139 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2140 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2141 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2142 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2144 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2147 ! For the GB potential convert sigma'**2 into chi'
2150 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2154 write (iout,'(/a/)') 'Parameters of the BP potential:'
2155 write (iout,'(a/)') 'The epsilon array:'
2156 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2157 write (iout,'(/a)') 'One-body parameters:'
2158 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2160 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2161 chip(i),alp(i),i=1,ntyp)
2164 !--------------------- GBV potential -----------------------------------
2166 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2167 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2168 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2169 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2171 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2172 write (iout,'(a/)') 'The epsilon array:'
2173 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2174 write (iout,'(/a)') 'One-body parameters:'
2175 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2176 's||/s_|_^2',' chip ',' alph '
2177 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2178 sigii(i),chip(i),alp(i),i=1,ntyp)
2181 write(iout,*)"Wrong ipot"
2187 !-----------------------------------------------------------------------
2188 ! Calculate the "working" parameters of SC interactions.
2190 !el from module energy - COMMON.INTERACT-------
2191 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2192 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2193 allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2194 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2196 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2209 sc_aa_tube_par(:)=0.0d0
2210 sc_bb_tube_par(:)=0.0d0
2212 !--------------------------------
2217 epslip(i,j)=epslip(j,i)
2222 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2223 sigma(j,i)=sigma(i,j)
2224 rs0(i,j)=dwa16*sigma(i,j)
2228 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2229 'Working parameters of the SC interactions:',&
2230 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2235 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2244 sigeps=dsign(1.0D0,epsij)
2246 aa_aq(i,j)=epsij*rrij*rrij
2247 bb_aq(i,j)=-sigeps*epsij*rrij
2248 aa_aq(j,i)=aa_aq(i,j)
2249 bb_aq(j,i)=bb_aq(i,j)
2250 epsijlip=epslip(i,j)
2251 sigeps=dsign(1.0D0,epsijlip)
2252 epsijlip=dabs(epsijlip)
2253 aa_lip(i,j)=epsijlip*rrij*rrij
2254 bb_lip(i,j)=-sigeps*epsijlip*rrij
2255 aa_lip(j,i)=aa_lip(i,j)
2256 bb_lip(j,i)=bb_lip(i,j)
2257 !C write(iout,*) aa_lip
2259 sigt1sq=sigma0(i)**2
2260 sigt2sq=sigma0(j)**2
2263 ratsig1=sigt2sq/sigt1sq
2264 ratsig2=1.0D0/ratsig1
2265 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2266 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2267 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2271 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2272 sigmaii(i,j)=rsum_max
2273 sigmaii(j,i)=rsum_max
2275 ! sigmaii(i,j)=r0(i,j)
2276 ! sigmaii(j,i)=r0(i,j)
2278 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2279 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2280 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2281 augm(i,j)=epsij*r_augm**(2*expon)
2282 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2289 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2290 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2291 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2296 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2297 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2298 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2299 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2300 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2301 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2302 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2303 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2304 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2305 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2306 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2307 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2308 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2309 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2310 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2311 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2320 read (isidep_nucl,*) ipot_nucl
2321 ! print *,"TU?!",ipot_nucl
2322 if (ipot_nucl.eq.1) then
2323 do i=1,ntyp_molec(2)
2324 do j=i,ntyp_molec(2)
2325 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2326 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2330 do i=1,ntyp_molec(2)
2331 do j=i,ntyp_molec(2)
2332 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2333 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2334 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2338 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2339 do i=1,ntyp_molec(2)
2340 do j=i,ntyp_molec(2)
2341 rrij=sigma_nucl(i,j)
2345 epsij=4*eps_nucl(i,j)
2346 sigeps=dsign(1.0D0,epsij)
2348 aa_nucl(i,j)=epsij*rrij*rrij
2349 bb_nucl(i,j)=-sigeps*epsij*rrij
2350 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2351 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2352 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2353 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2354 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2355 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2358 aa_nucl(i,j)=aa_nucl(j,i)
2359 bb_nucl(i,j)=bb_nucl(j,i)
2360 ael3_nucl(i,j)=ael3_nucl(j,i)
2361 ael6_nucl(i,j)=ael6_nucl(j,i)
2362 ael63_nucl(i,j)=ael63_nucl(j,i)
2363 ael32_nucl(i,j)=ael32_nucl(j,i)
2364 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2365 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2366 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2367 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2368 eps_nucl(i,j)=eps_nucl(j,i)
2369 sigma_nucl(i,j)=sigma_nucl(j,i)
2370 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2374 write(iout,*) "tube param"
2375 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2376 ccavtubpep,dcavtubpep,tubetranenepep
2377 sigmapeptube=sigmapeptube**6
2378 sigeps=dsign(1.0D0,epspeptube)
2379 epspeptube=dabs(epspeptube)
2380 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2381 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2382 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2384 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2385 ccavtub(i),dcavtub(i),tubetranene(i)
2386 sigmasctube=sigmasctube**6
2387 sigeps=dsign(1.0D0,epssctube)
2388 epssctube=dabs(epssctube)
2389 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2390 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2391 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2394 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2399 ! Define the SC-p interaction constants (hard-coded; old style)
2402 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2404 ! aad(i,1)=0.3D0*4.0D0**12
2405 ! Following line for constants currently implemented
2406 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2407 aad(i,1)=1.5D0*4.0D0**12
2408 ! aad(i,1)=0.17D0*5.6D0**12
2410 ! "Soft" SC-p repulsion
2412 ! Following line for constants currently implemented
2413 ! aad(i,1)=0.3D0*4.0D0**6
2414 ! "Hard" SC-p repulsion
2415 bad(i,1)=3.0D0*4.0D0**6
2416 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2425 ! 8/9/01 Read the SC-p interaction constants from file
2428 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2431 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2432 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2433 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2434 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2438 write (iout,*) "Parameters of SC-p interactions:"
2440 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2441 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2446 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2448 do i=1,ntyp_molec(2)
2449 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
2451 do i=1,ntyp_molec(2)
2452 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2453 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2455 r0pp=1.12246204830937298142*5.16158
2461 ! Define the constants of the disulfide bridge
2465 ! Old arbitrary potential - commented out.
2470 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2471 ! energy surface of diethyl disulfide.
2472 ! A. Liwo and U. Kozlowska, 11/24/03
2489 write (iout,'(/a)') "Disulfide bridge parameters:"
2490 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2491 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2492 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2493 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2497 111 write (iout,*) "Error reading bending energy parameters."
2499 112 write (iout,*) "Error reading rotamer energy parameters."
2501 113 write (iout,*) "Error reading torsional energy parameters."
2503 114 write (iout,*) "Error reading double torsional energy parameters."
2505 115 write (iout,*) &
2506 "Error reading cumulant (multibody energy) parameters."
2508 116 write (iout,*) "Error reading electrostatic energy parameters."
2510 117 write (iout,*) "Error reading side chain interaction parameters."
2512 118 write (iout,*) "Error reading SCp interaction parameters."
2514 119 write (iout,*) "Error reading SCCOR parameters"
2517 call MPI_Finalize(Ierror)
2521 end subroutine parmread
2523 !-----------------------------------------------------------------------------
2525 !-----------------------------------------------------------------------------
2526 subroutine printmat(ldim,m,n,iout,key,a)
2529 character(len=3),dimension(n) :: key
2530 real(kind=8),dimension(ldim,n) :: a
2532 integer :: i,j,k,m,iout,nlim
2536 write (iout,1000) (key(k),k=i,nlim)
2538 1000 format (/5x,8(6x,a3))
2539 1020 format (/80(1h-)/)
2541 write (iout,1010) key(j),(a(j,k),k=i,nlim)
2544 1010 format (a3,2x,8(f9.4))
2546 end subroutine printmat
2547 !-----------------------------------------------------------------------------
2549 !-----------------------------------------------------------------------------
2551 ! Read the PDB file and convert the peptide geometry into virtual-chain
2554 use energy_data, only: itype,istype
2558 use control, only: rescode,sugarcode
2559 ! implicit real*8 (a-h,o-z)
2560 ! include 'DIMENSIONS'
2561 ! include 'COMMON.LOCAL'
2562 ! include 'COMMON.VAR'
2563 ! include 'COMMON.CHAIN'
2564 ! include 'COMMON.INTERACT'
2565 ! include 'COMMON.IOUNITS'
2566 ! include 'COMMON.GEO'
2567 ! include 'COMMON.NAMES'
2568 ! include 'COMMON.CONTROL'
2569 ! include 'COMMON.DISTFIT'
2570 ! include 'COMMON.SETUP'
2571 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
2573 logical :: lprn=.true.,fail
2574 real(kind=8),dimension(3) :: e1,e2,e3
2575 real(kind=8) :: dcj,efree_temp
2576 character(len=3) :: seq,res
2577 character(len=5) :: atom
2578 character(len=80) :: card
2579 real(kind=8),dimension(3,20) :: sccor
2580 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
2581 integer :: isugar,molecprev
2582 character*1 :: sugar
2584 real(kind=8),dimension(3) :: ccc
2586 integer,dimension(2,maxres/3) :: hfrag_alloc
2587 integer,dimension(4,maxres/3) :: bfrag_alloc
2588 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2589 real(kind=8),dimension(:,:), allocatable :: c_temporary
2590 integer,dimension(:,:) , allocatable :: itype_temporary
2591 integer,dimension(:),allocatable :: istype_temp
2598 ! write (2,*) "UNRES_PDB",unres_pdb
2606 !-----------------------------
2607 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2608 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2611 read (ipdbin,'(a80)',end=10) card
2612 ! write (iout,'(a)') card
2613 if (card(:5).eq.'HELIX') then
2616 read(card(22:25),*) hfrag(1,nhfrag)
2617 read(card(34:37),*) hfrag(2,nhfrag)
2619 if (card(:5).eq.'SHEET') then
2622 read(card(24:26),*) bfrag(1,nbfrag)
2623 read(card(35:37),*) bfrag(2,nbfrag)
2624 !rc----------------------------------------
2625 !rc to be corrected !!!
2626 bfrag(3,nbfrag)=bfrag(1,nbfrag)
2627 bfrag(4,nbfrag)=bfrag(2,nbfrag)
2628 !rc----------------------------------------
2630 if (card(:3).eq.'END') then
2632 else if (card(:3).eq.'TER') then
2637 itype(ires_old,molecule)=ntyp1_molec(molecule)
2638 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
2639 nres_molec(molecule)=nres_molec(molecule)+2
2641 ! write (iout,*) "Chain ended",ires,ishift,ires_old
2644 dc(j,ires)=sccor(j,iii)
2647 call sccenter(ires,iii,sccor)
2653 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2654 ! Fish out the ATOM cards.
2655 ! write(iout,*) 'card',card(1:20)
2656 if (index(card(1:4),'ATOM').gt.0) then
2657 read (card(12:16),*) atom
2658 ! write (iout,*) "! ",atom," !",ires
2659 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
2660 read (card(23:26),*) ires
2661 read (card(18:20),'(a3)') res
2662 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
2663 ! & " ires_old",ires_old
2664 ! write (iout,*) "ishift",ishift," ishift1",ishift1
2665 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2666 if (ires-ishift+ishift1.ne.ires_old) then
2667 ! Calculate the CM of the preceding residue.
2668 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2670 ! write (iout,*) "Calculating sidechain center iii",iii
2673 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2676 call sccenter(ires_old,iii,sccor)
2680 ! Start new residue.
2681 if (res.eq.'Cl-' .or. res.eq.'Na+') then
2684 else if (ibeg.eq.1) then
2685 write (iout,*) "BEG ires",ires
2687 if (res.ne.'GLY' .and. res.ne. 'ACE') then
2690 nres_molec(molecule)=nres_molec(molecule)+1
2692 ires=ires-ishift+ishift1
2694 ! write (iout,*) "ishift",ishift," ires",ires,&
2695 ! " ires_old",ires_old
2697 else if (ibeg.eq.2) then
2699 ishift=-ires_old+ires-1 !!!!!
2700 ishift1=ishift1-1 !!!!!
2701 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2702 ires=ires-ishift+ishift1
2703 ! print *,ires,ishift,ishift1
2707 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2708 ires=ires-ishift+ishift1
2711 ! print *,'atom',ires,atom
2712 if (res.eq.'ACE' .or. res.eq.'NHE') then
2715 if (atom.eq.'CA '.or.atom.eq.'N ') then
2717 itype(ires,molecule)=rescode(ires,res,0,molecule)
2718 ! nres_molec(molecule)=nres_molec(molecule)+1
2721 itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
2722 ! nres_molec(molecule)=nres_molec(molecule)+1
2723 read (card(19:19),'(a1)') sugar
2724 isugar=sugarcode(sugar,ires)
2725 ! if (ibeg.eq.1) then
2729 ! print *,"ires=",ires,istype(ires)
2735 ires=ires-ishift+ishift1
2737 ! write (iout,*) "ires_old",ires_old," ires",ires
2738 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2741 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
2742 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2743 res.eq.'NHE'.and.atom(:2).eq.'HN') then
2744 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2745 ! print *,ires,ishift,ishift1
2746 ! write (iout,*) "backbone ",atom
2748 write (iout,'(2i3,2x,a,3f8.3)') &
2749 ires,itype(ires,1),res,(c(j,ires),j=1,3)
2752 nres_molec(molecule)=nres_molec(molecule)+1
2754 sccor(j,iii)=c(j,ires)
2756 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
2757 atom.eq."C2'" .or. atom.eq."C3'" &
2758 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
2759 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
2760 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
2761 ! print *,ires,ishift,ishift1
2765 ! sccor(j,iii)=c(j,ires)
2768 c(j,ires)=c(j,ires)+ccc(j)/5.0
2770 print *,counter,molecule
2771 if (counter.eq.5) then
2773 nres_molec(molecule)=nres_molec(molecule)+1
2775 ! sccor(j,iii)=c(j,ires)
2779 ! print *, "ATOM",atom(1:3)
2780 else if (atom.eq."C5'") then
2781 read (card(19:19),'(a1)') sugar
2782 isugar=sugarcode(sugar,ires)
2787 ! print *,ires,istype(ires)
2790 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2793 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2795 ! write (*,*) card(23:27),ires,itype(ires,1)
2796 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2797 atom.ne.'N' .and. atom.ne.'C' .and. &
2798 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2799 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
2800 .and. atom.ne.'P '.and. &
2801 atom(1:1).ne.'H' .and. &
2802 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
2804 ! write (iout,*) "sidechain ",atom
2805 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2806 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
2807 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2809 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2812 else if ((ions).and.(card(1:6).eq.'HETATM')) then
2814 read (card(12:16),*) atom
2815 print *,"HETATOM", atom
2816 read (card(18:20),'(a3)') res
2817 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
2818 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
2819 .or.(atom(1:2).eq.'K ')) &
2822 if (molecule.ne.5) molecprev=molecule
2824 nres_molec(molecule)=nres_molec(molecule)+1
2825 itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
2826 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2830 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
2831 if (ires.eq.0) return
2832 ! Calculate dummy residue coordinates inside the "chain" of a multichain
2835 if (ires_old.ne.ires) nres_molec(molecule)=nres_molec(molecule)-2
2836 ! print *,'I have', nres_molec(:)
2838 do k=1,4 ! ions are without dummy
2839 if (nres_molec(k).eq.0) cycle
2841 ! write (iout,*) i,itype(i,1)
2842 ! if (itype(i,1).eq.ntyp1) then
2843 ! write (iout,*) "dummy",i,itype(i,1)
2845 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
2846 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
2850 if (itype(i,k).eq.ntyp1_molec(k)) then
2851 if (itype(i+1,k).eq.ntyp1_molec(k)) then
2852 if (itype(i+2,k).eq.0) then
2853 ! print *,"masz sieczke"
2855 if (itype(i+2,j).ne.0) then
2857 itype(i+1,j)=ntyp1_molec(j)
2858 nres_molec(k)=nres_molec(k)-1
2859 nres_molec(j)=nres_molec(j)+1
2865 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
2866 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
2867 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
2869 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2870 ! print *,i,'tu dochodze'
2871 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
2879 c(j,i)=c(j,i-1)-1.9d0*e2(j)
2883 dcj=(c(j,i-2)-c(j,i-3))/2.0
2884 if (dcj.eq.0) dcj=1.23591524223
2889 else !itype(i+1,1).eq.ntyp1
2891 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2892 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
2899 c(j,i)=c(j,i+1)-1.9d0*e2(j)
2903 dcj=(c(j,i+3)-c(j,i+2))/2.0
2904 if (dcj.eq.0) dcj=1.23591524223
2909 endif !itype(i+1,1).eq.ntyp1
2910 endif !itype.eq.ntyp1
2914 ! Calculate the CM of the last side chain.
2918 dc(j,ires)=sccor(j,iii)
2921 call sccenter(ires,iii,sccor)
2927 ! print *,"molecule",molecule
2928 if ((itype(nres,1).ne.10)) then
2930 if (molecule.eq.5) molecule=molecprev
2931 itype(nres,molecule)=ntyp1_molec(molecule)
2932 nres_molec(molecule)=nres_molec(molecule)+1
2934 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2935 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
2942 c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
2946 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
2947 c(j,nres)=c(j,nres-1)+dcj
2948 c(j,2*nres)=c(j,nres)
2952 ! print *,'I have',nres, nres_molec(:)
2954 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
2956 if (nres.ne.nres0) then
2957 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
2959 stop "Error nres value in WHAM input"
2962 !---------------------------------
2963 !el reallocate tables
2966 ! hfrag_alloc(j,i)=hfrag(j,i)
2969 ! bfrag_alloc(j,i)=bfrag(j,i)
2975 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
2976 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
2977 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
2978 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
2982 ! hfrag(j,i)=hfrag_alloc(j,i)
2987 ! bfrag(j,i)=bfrag_alloc(j,i)
2990 !el end reallocate tables
2991 !---------------------------------
2999 c(j,2*nres)=c(j,nres)
3002 if (itype(1,1).eq.ntyp1) then
3006 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3007 call refsys(2,3,4,e1,e2,e3,fail)
3014 c(j,1)=c(j,2)-1.9d0*e2(j)
3018 dcj=(c(j,4)-c(j,3))/2.0
3024 ! First lets assign correct dummy to correct type of chain
3026 if (itype(1,1).eq.ntyp1) then
3027 if (itype(2,1).eq.0) then
3029 if (itype(2,j).ne.0) then
3031 itype(1,j)=ntyp1_molec(j)
3032 nres_molec(1)=nres_molec(1)-1
3033 nres_molec(j)=nres_molec(j)+1
3040 print *,'I have',nres, nres_molec(:)
3042 ! Copy the coordinates to reference coordinates
3048 ! Calculate internal coordinates.
3050 write (iout,'(/a)') &
3051 "Cartesian coordinates of the reference structure"
3052 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3053 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3055 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3056 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3057 (c(j,ires+nres),j=1,3)
3060 ! znamy już nres wiec mozna alokowac tablice
3061 ! Calculate internal coordinates.
3062 if(me.eq.king.or..not.out1file)then
3063 write (iout,'(a)') &
3064 "Backbone and SC coordinates as read from the PDB"
3066 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
3067 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3068 (c(j,nres+ires),j=1,3)
3071 ! NOW LETS ROCK! SORTING
3072 allocate(c_temporary(3,2*nres))
3073 allocate(itype_temporary(nres,5))
3074 allocate(molnum(nres))
3075 allocate(istype_temp(nres))
3076 itype_temporary(:,:)=0
3080 if (itype(i,k).ne.0) then
3082 c_temporary(j,seqalingbegin)=c(j,i)
3083 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3086 itype_temporary(seqalingbegin,k)=itype(i,k)
3087 istype_temp(seqalingbegin)=istype(i)
3088 molnum(seqalingbegin)=k
3089 seqalingbegin=seqalingbegin+1
3095 c(j,i)=c_temporary(j,i)
3100 itype(i,k)=itype_temporary(i,k)
3101 istype(i)=istype_temp(i)
3104 if (itype(1,1).eq.ntyp1) then
3108 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3109 call refsys(2,3,4,e1,e2,e3,fail)
3116 c(j,1)=c(j,2)-1.9d0*e2(j)
3120 dcj=(c(j,4)-c(j,3))/2.0
3128 write (iout,'(/a)') &
3129 "Cartesian coordinates of the reference structure after sorting"
3130 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3131 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3133 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3134 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3135 (c(j,ires+nres),j=1,3)
3139 ! print *,seqalingbegin,nres
3140 if(.not.allocated(vbld)) then
3141 allocate(vbld(2*nres))
3146 if(.not.allocated(vbld_inv)) then
3147 allocate(vbld_inv(2*nres))
3153 if(.not.allocated(theta)) then
3154 allocate(theta(nres+2))
3158 if(.not.allocated(phi)) allocate(phi(nres+2))
3159 if(.not.allocated(alph)) allocate(alph(nres+2))
3160 if(.not.allocated(omeg)) allocate(omeg(nres+2))
3161 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3162 if(.not.allocated(phiref)) allocate(phiref(nres+2))
3163 if(.not.allocated(costtab)) allocate(costtab(nres))
3164 if(.not.allocated(sinttab)) allocate(sinttab(nres))
3165 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3166 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3167 if(.not.allocated(xxref)) allocate(xxref(nres))
3168 if(.not.allocated(yyref)) allocate(yyref(nres))
3169 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3170 if(.not.allocated(dc_norm)) then
3171 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3172 allocate(dc_norm(3,0:2*nres+2))
3176 call int_from_cart(.true.,.false.)
3177 call sc_loc_geom(.false.)
3179 thetaref(i)=theta(i)
3189 dc(j,i)=c(j,i+1)-c(j,i)
3190 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3195 dc(j,i+nres)=c(j,i+nres)-c(j,i)
3196 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3198 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3202 ! Copy the coordinates to reference coordinates
3203 ! Splits to single chain if occurs
3207 ! cref(j,i,cou)=c(j,i)
3211 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
3212 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
3213 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3214 !-----------------------------
3218 write (iout,*) "symetr", symetr
3221 !c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3223 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
3226 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3231 cref(j,i,cou)=c(j,i)
3232 cref(j,i+nres,cou)=c(j,i+nres)
3234 chain_rep(j,lll,kkk)=c(j,i)
3235 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
3239 write (iout,*) chain_length
3240 if (chain_length.eq.0) chain_length=nres
3242 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
3243 chain_rep(j,chain_length+nres,symetr) &
3244 =chain_rep(j,chain_length+nres,1)
3247 ! write (iout,*) "spraw lancuchy",chain_length,symetr
3249 ! do kkk=1,chain_length
3250 ! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
3254 ! makes copy of chains
3255 write (iout,*) "symetr", symetr
3260 if (symetr.gt.1) then
3267 write(iout,*) (tabperm(i,kkk),kkk=1,4)
3273 ! write (iout,*) i,icha
3274 do lll=1,chain_length
3276 if (cou.le.nres) then
3278 kupa=mod(lll,chain_length)
3279 iprzes=(kkk-1)*chain_length+lll
3280 if (kupa.eq.0) kupa=chain_length
3281 ! write (iout,*) "kupa", kupa
3282 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
3283 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
3290 !-koniec robienia kopii
3293 write (iout,*) "nowa struktura", nperm
3295 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
3297 cref(3,i,kkk),cref(1,nres+i,kkk),&
3298 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
3300 100 format (//' alpha-carbon coordinates ',&
3301 ' centroid coordinates'/ &
3302 ' ', 6X,'X',11X,'Y',11X,'Z', &
3303 10X,'X',11X,'Y',11X,'Z')
3304 110 format (a,'(',i3,')',6f12.5)
3310 bfrag(i,j)=bfrag(i,j)-ishift
3316 hfrag(i,j)=hfrag(i,j)-ishift
3322 end subroutine readpdb
3323 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3324 !-----------------------------------------------------------------------------
3326 !-----------------------------------------------------------------------------
3327 subroutine read_control
3341 use random, only: random_init
3342 ! implicit real*8 (a-h,o-z)
3343 ! include 'DIMENSIONS'
3345 use prng, only:prng_restart
3347 logical :: OKRandom!, prng_restart
3350 ! include 'COMMON.IOUNITS'
3351 ! include 'COMMON.TIME1'
3352 ! include 'COMMON.THREAD'
3353 ! include 'COMMON.SBRIDGE'
3354 ! include 'COMMON.CONTROL'
3355 ! include 'COMMON.MCM'
3356 ! include 'COMMON.MAP'
3357 ! include 'COMMON.HEADER'
3358 ! include 'COMMON.CSA'
3359 ! include 'COMMON.CHAIN'
3360 ! include 'COMMON.MUCA'
3361 ! include 'COMMON.MD'
3362 ! include 'COMMON.FFIELD'
3363 ! include 'COMMON.INTERACT'
3364 ! include 'COMMON.SETUP'
3365 !el integer :: KDIAG,ICORFL,IXDR
3366 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
3367 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
3368 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
3369 ! character(len=80) :: ucase
3370 character(len=640) :: controlcard
3372 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
3378 read (INP,'(a)') titel
3379 call card_concat(controlcard,.true.)
3380 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
3381 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
3382 call reada(controlcard,'SEED',seed,0.0D0)
3383 call random_init(seed)
3384 ! Set up the time limit (caution! The time must be input in minutes!)
3385 read_cart=index(controlcard,'READ_CART').gt.0
3386 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3387 call readi(controlcard,'SYM',symetr,1)
3388 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
3389 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
3390 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
3391 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
3392 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
3393 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
3394 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
3395 call reada(controlcard,'DRMS',drms,0.1D0)
3396 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3397 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
3398 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
3399 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
3400 write (iout,'(a,f10.1)')'DRMS = ',drms
3401 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
3402 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
3404 call readi(controlcard,'NZ_START',nz_start,0)
3405 call readi(controlcard,'NZ_END',nz_end,0)
3406 call readi(controlcard,'IZ_SC',iz_sc,0)
3407 timlim=60.0D0*timlim
3408 safety = 60.0d0*safety
3411 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3412 !C SHIELD keyword sets if the shielding effect of side-chains is used
3413 !C 0 denots no shielding is used all peptide are equally despite the
3414 !C solvent accesible area
3415 !C 1 the newly introduced function
3416 !C 2 reseved for further possible developement
3417 call readi(controlcard,'SHIELD',shield_mode,0)
3418 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3419 write(iout,*) "shield_mode",shield_mode
3420 !C Varibles set size of box
3421 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
3422 protein=index(controlcard,"PROTEIN").gt.0
3423 ions=index(controlcard,"IONS").gt.0
3424 nucleic=index(controlcard,"NUCLEIC").gt.0
3425 write (iout,*) "with_theta_constr ",with_theta_constr
3426 AFMlog=(index(controlcard,'AFM'))
3427 selfguide=(index(controlcard,'SELFGUIDE'))
3428 print *,'AFMlog',AFMlog,selfguide,"KUPA"
3429 call readi(controlcard,'GENCONSTR',genconstr,0)
3430 call reada(controlcard,'BOXX',boxxsize,100.0d0)
3431 call reada(controlcard,'BOXY',boxysize,100.0d0)
3432 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3433 call readi(controlcard,'TUBEMOD',tubemode,0)
3434 write (iout,*) TUBEmode,"TUBEMODE"
3435 if (TUBEmode.gt.0) then
3436 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3437 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3438 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3439 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3440 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3441 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3442 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3443 buftubebot=bordtubebot+tubebufthick
3444 buftubetop=bordtubetop-tubebufthick
3447 ! CUTOFFF ON ELECTROSTATICS
3448 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
3449 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3450 write(iout,*) "R_CUT_ELE=",r_cut_ele
3451 ! Lipidic parameters
3452 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3453 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3454 if (lipthick.gt.0.0d0) then
3455 bordliptop=(boxzsize+lipthick)/2.0
3456 bordlipbot=bordliptop-lipthick
3457 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3458 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3459 buflipbot=bordlipbot+lipbufthick
3460 bufliptop=bordliptop-lipbufthick
3461 if ((lipbufthick*2.0d0).gt.lipthick) &
3462 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3463 endif !lipthick.gt.0
3464 write(iout,*) "bordliptop=",bordliptop
3465 write(iout,*) "bordlipbot=",bordlipbot
3466 write(iout,*) "bufliptop=",bufliptop
3467 write(iout,*) "buflipbot=",buflipbot
3468 write (iout,*) "SHIELD MODE",shield_mode
3470 !C-------------------------
3471 minim=(index(controlcard,'MINIMIZE').gt.0)
3472 dccart=(index(controlcard,'CART').gt.0)
3473 overlapsc=(index(controlcard,'OVERLAP').gt.0)
3474 overlapsc=.not.overlapsc
3475 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3476 searchsc=.not.searchsc
3477 sideadd=(index(controlcard,'SIDEADD').gt.0)
3478 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3479 outpdb=(index(controlcard,'PDBOUT').gt.0)
3480 outmol2=(index(controlcard,'MOL2OUT').gt.0)
3481 pdbref=(index(controlcard,'PDBREF').gt.0)
3482 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3483 indpdb=index(controlcard,'PDBSTART')
3484 extconf=(index(controlcard,'EXTCONF').gt.0)
3485 call readi(controlcard,'IPRINT',iprint,0)
3486 call readi(controlcard,'MAXGEN',maxgen,10000)
3487 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3488 call readi(controlcard,"KDIAG",kdiag,0)
3489 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3490 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3491 write (iout,*) "RESCALE_MODE",rescale_mode
3492 split_ene=index(controlcard,'SPLIT_ENE').gt.0
3493 if (index(controlcard,'REGULAR').gt.0.0D0) then
3494 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3498 if (index(controlcard,'CHECKGRAD').gt.0) then
3500 if (index(controlcard,'CART').gt.0) then
3502 elseif (index(controlcard,'CARINT').gt.0) then
3507 elseif (index(controlcard,'THREAD').gt.0) then
3509 call readi(controlcard,'THREAD',nthread,0)
3510 if (nthread.gt.0) then
3511 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3514 write (iout,'(a)')'A number has to follow the THREAD keyword.'
3515 stop 'Error termination in Read_Control.'
3517 else if (index(controlcard,'MCMA').gt.0) then
3519 else if (index(controlcard,'MCEE').gt.0) then
3521 else if (index(controlcard,'MULTCONF').gt.0) then
3523 else if (index(controlcard,'MAP').gt.0) then
3525 call readi(controlcard,'MAP',nmap,0)
3526 else if (index(controlcard,'CSA').gt.0) then
3528 !rc else if (index(controlcard,'ZSCORE').gt.0) then
3530 !rc ZSCORE is rm from UNRES, modecalc=9 is available
3533 !fcm else if (index(controlcard,'MCMF').gt.0) then
3535 else if (index(controlcard,'SOFTREG').gt.0) then
3537 else if (index(controlcard,'CHECK_BOND').gt.0) then
3539 else if (index(controlcard,'TEST').gt.0) then
3541 else if (index(controlcard,'MD').gt.0) then
3543 else if (index(controlcard,'RE ').gt.0) then
3547 lmuca=index(controlcard,'MUCA').gt.0
3548 call readi(controlcard,'MUCADYN',mucadyn,0)
3549 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3550 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3552 write (iout,*) 'MUCADYN=',mucadyn
3553 write (iout,*) 'MUCASMOOTH=',muca_smooth
3556 iscode=index(controlcard,'ONE_LETTER')
3557 indphi=index(controlcard,'PHI')
3558 indback=index(controlcard,'BACK')
3559 iranconf=index(controlcard,'RAND_CONF')
3560 i2ndstr=index(controlcard,'USE_SEC_PRED')
3561 gradout=index(controlcard,'GRADOUT').gt.0
3562 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3563 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3564 if (me.eq.king .or. .not.out1file ) &
3565 write (iout,*) "DISTCHAINMAX",distchainmax
3567 if(me.eq.king.or..not.out1file) &
3568 write (iout,'(2a)') diagmeth(kdiag),&
3569 ' routine used to diagonalize matrices.'
3570 if (shield_mode.gt.0) then
3572 !C VSolvSphere the volume of solving sphere
3574 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3575 !C there will be no distinction between proline peptide group and normal peptide
3576 !C group in case of shielding parameters
3577 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3578 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3579 write (iout,*) VSolvSphere,VSolvSphere_div
3580 !C long axis of side chain
3582 long_r_sidechain(i)=vbldsc0(1,i)
3583 short_r_sidechain(i)=sigma0(i)
3584 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3590 end subroutine read_control
3591 !-----------------------------------------------------------------------------
3592 subroutine read_REMDpar
3594 ! Read REMD settings
3601 use control_data, only:out1file
3602 ! implicit real*8 (a-h,o-z)
3603 ! include 'DIMENSIONS'
3604 ! include 'COMMON.IOUNITS'
3605 ! include 'COMMON.TIME1'
3606 ! include 'COMMON.MD'
3609 !el include 'COMMON.LANGEVIN'
3611 !el include 'COMMON.LANGEVIN.lang0'
3613 ! include 'COMMON.INTERACT'
3614 ! include 'COMMON.NAMES'
3615 ! include 'COMMON.GEO'
3616 ! include 'COMMON.REMD'
3617 ! include 'COMMON.CONTROL'
3618 ! include 'COMMON.SETUP'
3619 ! character(len=80) :: ucase
3620 character(len=320) :: controlcard
3621 character(len=3200) :: controlcard1
3622 integer :: iremd_m_total
3625 ! real(kind=8) :: var,ene
3627 if(me.eq.king.or..not.out1file) &
3628 write (iout,*) "REMD setup"
3630 call card_concat(controlcard,.true.)
3631 call readi(controlcard,"NREP",nrep,3)
3632 call readi(controlcard,"NSTEX",nstex,1000)
3633 call reada(controlcard,"RETMIN",retmin,10.0d0)
3634 call reada(controlcard,"RETMAX",retmax,1000.0d0)
3635 mremdsync=(index(controlcard,'SYNC').gt.0)
3636 call readi(controlcard,"NSYN",i_sync_step,100)
3637 restart1file=(index(controlcard,'REST1FILE').gt.0)
3638 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3639 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3640 if(max_cache_traj_use.gt.max_cache_traj) &
3641 max_cache_traj_use=max_cache_traj
3642 if(me.eq.king.or..not.out1file) then
3643 !d if (traj1file) then
3644 !rc caching is in testing - NTWX is not ignored
3645 !d write (iout,*) "NTWX value is ignored"
3646 !d write (iout,*) " trajectory is stored to one file by master"
3647 !d write (iout,*) " before exchange at NSTEX intervals"
3649 write (iout,*) "NREP= ",nrep
3650 write (iout,*) "NSTEX= ",nstex
3651 write (iout,*) "SYNC= ",mremdsync
3652 write (iout,*) "NSYN= ",i_sync_step
3653 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3656 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3657 if (index(controlcard,'TLIST').gt.0) then
3659 call card_concat(controlcard1,.true.)
3660 read(controlcard1,*) (remd_t(i),i=1,nrep)
3661 if(me.eq.king.or..not.out1file) &
3662 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3665 if (index(controlcard,'MLIST').gt.0) then
3667 call card_concat(controlcard1,.true.)
3668 read(controlcard1,*) (remd_m(i),i=1,nrep)
3669 if(me.eq.king.or..not.out1file) then
3670 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3673 iremd_m_total=iremd_m_total+remd_m(i)
3675 write (iout,*) 'Total number of replicas ',iremd_m_total
3678 if(me.eq.king.or..not.out1file) &
3679 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3681 end subroutine read_REMDpar
3682 !-----------------------------------------------------------------------------
3683 subroutine read_MDpar
3687 use control_data, only: r_cut,rlamb,out1file
3689 use geometry_data, only: pi
3691 ! implicit real*8 (a-h,o-z)
3692 ! include 'DIMENSIONS'
3693 ! include 'COMMON.IOUNITS'
3694 ! include 'COMMON.TIME1'
3695 ! include 'COMMON.MD'
3698 !el include 'COMMON.LANGEVIN'
3700 !el include 'COMMON.LANGEVIN.lang0'
3702 ! include 'COMMON.INTERACT'
3703 ! include 'COMMON.NAMES'
3704 ! include 'COMMON.GEO'
3705 ! include 'COMMON.SETUP'
3706 ! include 'COMMON.CONTROL'
3707 ! include 'COMMON.SPLITELE'
3708 ! character(len=80) :: ucase
3709 character(len=320) :: controlcard
3714 call card_concat(controlcard,.true.)
3715 call readi(controlcard,"NSTEP",n_timestep,1000000)
3716 call readi(controlcard,"NTWE",ntwe,100)
3717 call readi(controlcard,"NTWX",ntwx,1000)
3718 call reada(controlcard,"DT",d_time,1.0d-1)
3719 call reada(controlcard,"DVMAX",dvmax,2.0d1)
3720 call reada(controlcard,"DAMAX",damax,1.0d1)
3721 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3722 call readi(controlcard,"LANG",lang,0)
3723 RESPA = index(controlcard,"RESPA") .gt. 0
3724 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3725 ntime_split0=ntime_split
3726 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3727 ntime_split0=ntime_split
3728 call reada(controlcard,"R_CUT",r_cut,2.0d0)
3729 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3730 rest = index(controlcard,"REST").gt.0
3731 tbf = index(controlcard,"TBF").gt.0
3732 usampl = index(controlcard,"USAMPL").gt.0
3733 mdpdb = index(controlcard,"MDPDB").gt.0
3734 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3735 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3736 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3737 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3738 if (count_reset_moment.eq.0) count_reset_moment=1000000000
3739 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3740 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3741 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3742 if (count_reset_vel.eq.0) count_reset_vel=1000000000
3743 large = index(controlcard,"LARGE").gt.0
3744 print_compon = index(controlcard,"PRINT_COMPON").gt.0
3745 rattle = index(controlcard,"RATTLE").gt.0
3746 ! if performing umbrella sampling, fragments constrained are read from the fragment file
3752 if(me.eq.king.or..not.out1file) then
3754 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3756 write (iout,'(a)') "The units are:"
3757 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3758 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3759 " acceleration: angstrom/(48.9 fs)**2"
3760 write (iout,'(a)') "energy: kcal/mol, temperature: K"
3762 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3763 write (iout,'(a60,f10.5,a)') &
3764 "Initial time step of numerical integration:",d_time,&
3766 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3768 write (iout,'(2a,i4,a)') &
3769 "A-MTS algorithm used; initial time step for fast-varying",&
3770 " short-range forces split into",ntime_split," steps."
3771 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3772 r_cut," lambda",rlamb
3774 write (iout,'(2a,f10.5)') &
3775 "Maximum acceleration threshold to reduce the time step",&
3776 "/increase split number:",damax
3777 write (iout,'(2a,f10.5)') &
3778 "Maximum predicted energy drift to reduce the timestep",&
3779 "/increase split number:",edriftmax
3780 write (iout,'(a60,f10.5)') &
3781 "Maximum velocity threshold to reduce velocities:",dvmax
3782 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3783 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3784 if (rattle) write (iout,'(a60)') &
3785 "Rattle algorithm used to constrain the virtual bonds"
3789 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3790 call reada(controlcard,"RWAT",rwat,1.4d0)
3791 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3792 surfarea=index(controlcard,"SURFAREA").gt.0
3793 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3794 if(me.eq.king.or..not.out1file)then
3795 write (iout,'(/a,$)') "Langevin dynamics calculation"
3797 write (iout,'(a/)') &
3798 " with direct integration of Langevin equations"
3799 else if (lang.eq.2) then
3800 write (iout,'(a/)') " with TINKER stochasic MD integrator"
3801 else if (lang.eq.3) then
3802 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3803 else if (lang.eq.4) then
3804 write (iout,'(a/)') " in overdamped mode"
3806 write (iout,'(//a,i5)') &
3807 "=========== ERROR: Unknown Langevin dynamics mode:",lang
3810 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3811 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3812 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3813 write (iout,'(a60,f10.5)') &
3814 "Scaling factor of the friction forces:",scal_fric
3815 if (surfarea) write (iout,'(2a,i10,a)') &
3816 "Friction coefficients will be scaled by solvent-accessible",&
3817 " surface area every",reset_fricmat," steps."
3819 ! Calculate friction coefficients and bounds of stochastic forces
3820 eta=6*pi*cPoise*etawat
3821 if(me.eq.king.or..not.out1file) &
3822 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3825 do j=1,5 !types of molecules
3826 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
3827 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
3829 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
3830 do j=1,5 !types of molecules
3832 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
3833 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
3837 if(me.eq.king.or..not.out1file)then
3838 write (iout,'(/2a/)') &
3839 "Radii of site types and friction coefficients and std's of",&
3840 " stochastic forces of fully exposed sites"
3841 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
3843 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
3844 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
3848 if(me.eq.king.or..not.out1file)then
3849 write (iout,'(a)') "Berendsen bath calculation"
3850 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3851 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3853 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3854 count_reset_moment," steps"
3856 write (iout,'(a,i10,a)') &
3857 "Velocities will be reset at random every",count_reset_vel,&
3861 if(me.eq.king.or..not.out1file) &
3862 write (iout,'(a31)') "Microcanonical mode calculation"
3864 if(me.eq.king.or..not.out1file)then
3865 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3867 write(iout,*) "MD running with constraints."
3868 write(iout,*) "Equilibration time ", eq_time, " mtus."
3869 write(iout,*) "Constraining ", nfrag," fragments."
3870 write(iout,*) "Length of each fragment, weight and q0:"
3872 write (iout,*) "Set of restraints #",iset
3874 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3875 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3877 write(iout,*) "constraints between ", npair, "fragments."
3878 write(iout,*) "constraint pairs, weights and q0:"
3880 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3881 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3883 write(iout,*) "angle constraints within ", nfrag_back,&
3884 "backbone fragments."
3885 write(iout,*) "fragment, weights:"
3887 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3888 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3889 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3892 iset=mod(kolor,nset)+1
3895 if(me.eq.king.or..not.out1file) &
3896 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3898 end subroutine read_MDpar
3899 !-----------------------------------------------------------------------------
3903 ! implicit real*8 (a-h,o-z)
3904 ! include 'DIMENSIONS'
3905 ! include 'COMMON.MAP'
3906 ! include 'COMMON.IOUNITS'
3907 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3908 character(len=80) :: mapcard !,ucase
3911 ! real(kind=8) :: var,ene
3914 read (inp,'(a)') mapcard
3915 mapcard=ucase(mapcard)
3916 if (index(mapcard,'PHI').gt.0) then
3918 else if (index(mapcard,'THE').gt.0) then
3920 else if (index(mapcard,'ALP').gt.0) then
3922 else if (index(mapcard,'OME').gt.0) then
3925 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
3926 stop 'Error - illegal variable spec in MAP card.'
3928 call readi (mapcard,'RES1',res1(imap),0)
3929 call readi (mapcard,'RES2',res2(imap),0)
3930 if (res1(imap).eq.0) then
3931 res1(imap)=res2(imap)
3932 else if (res2(imap).eq.0) then
3933 res2(imap)=res1(imap)
3935 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
3936 write (iout,'(a)') &
3937 'Error - illegal definition of variable group in MAP.'
3938 stop 'Error - illegal definition of variable group in MAP.'
3940 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
3941 call reada(mapcard,'TO',ang_to(imap),0.0D0)
3942 call readi(mapcard,'NSTEP',nstep(imap),0)
3943 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
3944 write (iout,'(a)') &
3945 'Illegal boundary and/or step size specification in MAP.'
3946 stop 'Illegal boundary and/or step size specification in MAP.'
3950 end subroutine map_read
3951 !-----------------------------------------------------------------------------
3954 use control_data, only: vdisulf
3956 ! implicit real*8 (a-h,o-z)
3957 ! include 'DIMENSIONS'
3958 ! include 'COMMON.IOUNITS'
3959 ! include 'COMMON.GEO'
3960 ! include 'COMMON.CSA'
3961 ! include 'COMMON.BANK'
3962 ! include 'COMMON.CONTROL'
3963 ! character(len=80) :: ucase
3964 character(len=620) :: mcmcard
3966 ! integer :: ntf,ik,iw_pdb
3967 ! real(kind=8) :: var,ene
3969 call card_concat(mcmcard,.true.)
3971 call readi(mcmcard,'NCONF',nconf,50)
3972 call readi(mcmcard,'NADD',nadd,0)
3973 call readi(mcmcard,'JSTART',jstart,1)
3974 call readi(mcmcard,'JEND',jend,1)
3975 call readi(mcmcard,'NSTMAX',nstmax,500000)
3976 call readi(mcmcard,'N0',n0,1)
3977 call readi(mcmcard,'N1',n1,6)
3978 call readi(mcmcard,'N2',n2,4)
3979 call readi(mcmcard,'N3',n3,0)
3980 call readi(mcmcard,'N4',n4,0)
3981 call readi(mcmcard,'N5',n5,0)
3982 call readi(mcmcard,'N6',n6,10)
3983 call readi(mcmcard,'N7',n7,0)
3984 call readi(mcmcard,'N8',n8,0)
3985 call readi(mcmcard,'N9',n9,0)
3986 call readi(mcmcard,'N14',n14,0)
3987 call readi(mcmcard,'N15',n15,0)
3988 call readi(mcmcard,'N16',n16,0)
3989 call readi(mcmcard,'N17',n17,0)
3990 call readi(mcmcard,'N18',n18,0)
3992 vdisulf=(index(mcmcard,'DYNSS').gt.0)
3994 call readi(mcmcard,'NDIFF',ndiff,2)
3995 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
3996 call readi(mcmcard,'IS1',is1,1)
3997 call readi(mcmcard,'IS2',is2,8)
3998 call readi(mcmcard,'NRAN0',nran0,4)
3999 call readi(mcmcard,'NRAN1',nran1,2)
4000 call readi(mcmcard,'IRR',irr,1)
4001 call readi(mcmcard,'NSEED',nseed,20)
4002 call readi(mcmcard,'NTOTAL',ntotal,10000)
4003 call reada(mcmcard,'CUT1',cut1,2.0d0)
4004 call reada(mcmcard,'CUT2',cut2,5.0d0)
4005 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4006 call readi(mcmcard,'ICMAX',icmax,3)
4007 call readi(mcmcard,'IRESTART',irestart,0)
4008 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
4011 call reada(mcmcard,'DELE',dele,20.0d0)
4012 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4013 call readi(mcmcard,'IREF',iref,0)
4014 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4015 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4016 call readi(mcmcard,'NCONF_IN',nconf_in,0)
4017 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4018 write (iout,*) "NCONF_IN",nconf_in
4020 end subroutine csaread
4021 !-----------------------------------------------------------------------------
4025 use control_data, only: MaxMoveType
4028 ! implicit real*8 (a-h,o-z)
4029 ! include 'DIMENSIONS'
4030 ! include 'COMMON.MCM'
4031 ! include 'COMMON.MCE'
4032 ! include 'COMMON.IOUNITS'
4033 ! character(len=80) :: ucase
4034 character(len=320) :: mcmcard
4037 ! real(kind=8) :: var,ene
4039 call card_concat(mcmcard,.true.)
4040 call readi(mcmcard,'MAXACC',maxacc,100)
4041 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4042 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4043 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4044 call readi(mcmcard,'MAXREPM',maxrepm,200)
4045 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4046 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4047 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4048 call reada(mcmcard,'E_UP',e_up,5.0D0)
4049 call reada(mcmcard,'DELTE',delte,0.1D0)
4050 call readi(mcmcard,'NSWEEP',nsweep,5)
4051 call readi(mcmcard,'NSTEPH',nsteph,0)
4052 call readi(mcmcard,'NSTEPC',nstepc,0)
4053 call reada(mcmcard,'TMIN',tmin,298.0D0)
4054 call reada(mcmcard,'TMAX',tmax,298.0D0)
4055 call readi(mcmcard,'NWINDOW',nwindow,0)
4056 call readi(mcmcard,'PRINT_MC',print_mc,0)
4057 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4058 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4059 ent_read=(index(mcmcard,'ENT_READ').gt.0)
4060 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4061 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4062 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4063 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4064 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4065 if (nwindow.gt.0) then
4066 allocate(winstart(nwindow)) !!el (maxres)
4067 allocate(winend(nwindow)) !!el
4068 allocate(winlen(nwindow)) !!el
4069 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4071 winlen(i)=winend(i)-winstart(i)+1
4074 if (tmax.lt.tmin) tmax=tmin
4075 if (tmax.eq.tmin) then
4079 if (nstepc.gt.0 .and. nsteph.gt.0) then
4080 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4081 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4083 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4084 ! Probabilities of different move types
4085 sumpro_type(0)=0.0D0
4086 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4087 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4088 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4089 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
4090 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4091 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4092 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4094 print *,'i',i,' sumprotype',sumpro_type(i)
4095 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4096 print *,'i',i,' sumprotype',sumpro_type(i)
4099 end subroutine mcmread
4100 !-----------------------------------------------------------------------------
4101 subroutine read_minim
4104 ! implicit real*8 (a-h,o-z)
4105 ! include 'DIMENSIONS'
4106 ! include 'COMMON.MINIM'
4107 ! include 'COMMON.IOUNITS'
4108 ! character(len=80) :: ucase
4109 character(len=320) :: minimcard
4111 ! integer :: ntf,ik,iw_pdb
4112 ! real(kind=8) :: var,ene
4114 call card_concat(minimcard,.true.)
4115 call readi(minimcard,'MAXMIN',maxmin,2000)
4116 call readi(minimcard,'MAXFUN',maxfun,5000)
4117 call readi(minimcard,'MINMIN',minmin,maxmin)
4118 call readi(minimcard,'MINFUN',minfun,maxmin)
4119 call reada(minimcard,'TOLF',tolf,1.0D-2)
4120 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4121 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4122 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4123 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4124 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4125 'Options in energy minimization:'
4126 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4127 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4128 'MinMin:',MinMin,' MinFun:',MinFun,&
4129 ' TolF:',TolF,' RTolF:',RTolF
4131 end subroutine read_minim
4132 !-----------------------------------------------------------------------------
4133 subroutine openunits
4135 use MD_data, only: usampl
4138 use control_data, only:out1file
4139 use control, only: getenv_loc
4140 ! implicit real*8 (a-h,o-z)
4141 ! include 'DIMENSIONS'
4144 character(len=16) :: form,nodename
4145 integer :: nodelen,ierror,npos
4147 ! include 'COMMON.SETUP'
4148 ! include 'COMMON.IOUNITS'
4149 ! include 'COMMON.MD'
4150 ! include 'COMMON.CONTROL'
4151 integer :: lenpre,lenpot,lentmp !,ilen
4153 character(len=3) :: out1file_text !,ucase
4154 character(len=3) :: ll
4157 ! integer :: ntf,ik,iw_pdb
4158 ! real(kind=8) :: var,ene
4160 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4161 call getenv_loc("PREFIX",prefix)
4163 call getenv_loc("POT",pot)
4164 call getenv_loc("DIRTMP",tmpdir)
4165 call getenv_loc("CURDIR",curdir)
4166 call getenv_loc("OUT1FILE",out1file_text)
4167 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4168 out1file_text=ucase(out1file_text)
4169 if (out1file_text(1:1).eq."Y") then
4172 out1file=fg_rank.gt.0
4177 if (lentmp.gt.0) then
4178 write (*,'(80(1h!))')
4179 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
4180 write (*,'(80(1h!))')
4181 write (*,*)"All output files will be on node /tmp directory."
4183 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4184 if (me.eq.king) then
4185 write (*,*) "The master node is ",nodename
4186 else if (fg_rank.eq.0) then
4187 write (*,*) "I am the CG slave node ",nodename
4189 write (*,*) "I am the FG slave node ",nodename
4192 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
4193 lenpre = lentmp+lenpre+1
4195 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
4196 ! Get the names and open the input files
4197 #if defined(WINIFL) || defined(WINPGI)
4198 open(1,file=pref_orig(:ilen(pref_orig))// &
4199 '.inp',status='old',readonly,shared)
4200 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4201 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4202 ! Get parameter filenames and open the parameter files.
4203 call getenv_loc('BONDPAR',bondname)
4204 open (ibond,file=bondname,status='old',readonly,shared)
4205 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4206 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
4207 call getenv_loc('THETPAR',thetname)
4208 open (ithep,file=thetname,status='old',readonly,shared)
4209 call getenv_loc('ROTPAR',rotname)
4210 open (irotam,file=rotname,status='old',readonly,shared)
4211 call getenv_loc('TORPAR',torname)
4212 open (itorp,file=torname,status='old',readonly,shared)
4213 call getenv_loc('TORDPAR',tordname)
4214 open (itordp,file=tordname,status='old',readonly,shared)
4215 call getenv_loc('FOURIER',fouriername)
4216 open (ifourier,file=fouriername,status='old',readonly,shared)
4217 call getenv_loc('ELEPAR',elename)
4218 open (ielep,file=elename,status='old',readonly,shared)
4219 call getenv_loc('SIDEPAR',sidename)
4220 open (isidep,file=sidename,status='old',readonly,shared)
4222 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4223 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
4224 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4225 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
4226 call getenv_loc('TORPAR_NUCL',torname_nucl)
4227 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
4228 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4229 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
4230 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4231 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
4234 #elif (defined CRAY) || (defined AIX)
4235 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4237 ! print *,"Processor",myrank," opened file 1"
4238 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4239 ! print *,"Processor",myrank," opened file 9"
4240 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4241 ! Get parameter filenames and open the parameter files.
4242 call getenv_loc('BONDPAR',bondname)
4243 open (ibond,file=bondname,status='old',action='read')
4244 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4245 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4247 ! print *,"Processor",myrank," opened file IBOND"
4248 call getenv_loc('THETPAR',thetname)
4249 open (ithep,file=thetname,status='old',action='read')
4250 ! print *,"Processor",myrank," opened file ITHEP"
4251 call getenv_loc('ROTPAR',rotname)
4252 open (irotam,file=rotname,status='old',action='read')
4253 ! print *,"Processor",myrank," opened file IROTAM"
4254 call getenv_loc('TORPAR',torname)
4255 open (itorp,file=torname,status='old',action='read')
4256 ! print *,"Processor",myrank," opened file ITORP"
4257 call getenv_loc('TORDPAR',tordname)
4258 open (itordp,file=tordname,status='old',action='read')
4259 ! print *,"Processor",myrank," opened file ITORDP"
4260 call getenv_loc('SCCORPAR',sccorname)
4261 open (isccor,file=sccorname,status='old',action='read')
4262 ! print *,"Processor",myrank," opened file ISCCOR"
4263 call getenv_loc('FOURIER',fouriername)
4264 open (ifourier,file=fouriername,status='old',action='read')
4265 ! print *,"Processor",myrank," opened file IFOURIER"
4266 call getenv_loc('ELEPAR',elename)
4267 open (ielep,file=elename,status='old',action='read')
4268 ! print *,"Processor",myrank," opened file IELEP"
4269 call getenv_loc('SIDEPAR',sidename)
4270 open (isidep,file=sidename,status='old',action='read')
4272 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4273 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4274 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4275 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4276 call getenv_loc('TORPAR_NUCL',torname_nucl)
4277 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4278 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4279 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4280 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4281 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4283 call getenv_loc('LIPTRANPAR',liptranname)
4284 open (iliptranpar,file=liptranname,status='old',action='read')
4285 call getenv_loc('TUBEPAR',tubename)
4286 open (itube,file=tubename,status='old',action='read')
4288 ! print *,"Processor",myrank," opened file ISIDEP"
4289 ! print *,"Processor",myrank," opened parameter files"
4291 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
4292 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4293 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4294 ! Get parameter filenames and open the parameter files.
4295 call getenv_loc('BONDPAR',bondname)
4296 open (ibond,file=bondname,status='old')
4297 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4298 open (ibond_nucl,file=bondname_nucl,status='old')
4300 call getenv_loc('THETPAR',thetname)
4301 open (ithep,file=thetname,status='old')
4302 call getenv_loc('ROTPAR',rotname)
4303 open (irotam,file=rotname,status='old')
4304 call getenv_loc('TORPAR',torname)
4305 open (itorp,file=torname,status='old')
4306 call getenv_loc('TORDPAR',tordname)
4307 open (itordp,file=tordname,status='old')
4308 call getenv_loc('SCCORPAR',sccorname)
4309 open (isccor,file=sccorname,status='old')
4310 call getenv_loc('FOURIER',fouriername)
4311 open (ifourier,file=fouriername,status='old')
4312 call getenv_loc('ELEPAR',elename)
4313 open (ielep,file=elename,status='old')
4314 call getenv_loc('SIDEPAR',sidename)
4315 open (isidep,file=sidename,status='old')
4317 open (ithep_nucl,file=thetname_nucl,status='old')
4318 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4319 open (irotam_nucl,file=rotname_nucl,status='old')
4320 call getenv_loc('TORPAR_NUCL',torname_nucl)
4321 open (itorp_nucl,file=torname_nucl,status='old')
4322 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4323 open (itordp_nucl,file=tordname_nucl,status='old')
4324 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4325 open (isidep_nucl,file=sidename_nucl,status='old')
4327 call getenv_loc('LIPTRANPAR',liptranname)
4328 open (iliptranpar,file=liptranname,status='old')
4329 call getenv_loc('TUBEPAR',tubename)
4330 open (itube,file=tubename,status='old')
4332 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4334 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4335 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4336 ! Get parameter filenames and open the parameter files.
4337 call getenv_loc('BONDPAR',bondname)
4338 open (ibond,file=bondname,status='old',action='read')
4339 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4340 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4341 call getenv_loc('THETPAR',thetname)
4342 open (ithep,file=thetname,status='old',action='read')
4343 call getenv_loc('ROTPAR',rotname)
4344 open (irotam,file=rotname,status='old',action='read')
4345 call getenv_loc('TORPAR',torname)
4346 open (itorp,file=torname,status='old',action='read')
4347 call getenv_loc('TORDPAR',tordname)
4348 open (itordp,file=tordname,status='old',action='read')
4349 call getenv_loc('SCCORPAR',sccorname)
4350 open (isccor,file=sccorname,status='old',action='read')
4352 call getenv_loc('THETPARPDB',thetname_pdb)
4353 print *,"thetname_pdb ",thetname_pdb
4354 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
4355 print *,ithep_pdb," opened"
4357 call getenv_loc('FOURIER',fouriername)
4358 open (ifourier,file=fouriername,status='old',readonly)
4359 call getenv_loc('ELEPAR',elename)
4360 open (ielep,file=elename,status='old',readonly)
4361 call getenv_loc('SIDEPAR',sidename)
4362 open (isidep,file=sidename,status='old',readonly)
4364 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4365 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4366 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4367 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4368 call getenv_loc('TORPAR_NUCL',torname_nucl)
4369 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4370 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4371 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4372 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4373 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4375 call getenv_loc('LIPTRANPAR',liptranname)
4376 open (iliptranpar,file=liptranname,status='old',action='read')
4377 call getenv_loc('TUBEPAR',tubename)
4378 open (itube,file=tubename,status='old',action='read')
4381 call getenv_loc('ROTPARPDB',rotname_pdb)
4382 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
4385 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
4386 #if defined(WINIFL) || defined(WINPGI)
4387 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
4388 #elif (defined CRAY) || (defined AIX)
4389 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4391 open (iscpp_nucl,file=scpname_nucl,status='old')
4393 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4398 ! 8/9/01 In the newest version SCp interaction constants are read from a file
4399 ! Use -DOLDSCP to use hard-coded constants instead.
4401 call getenv_loc('SCPPAR',scpname)
4402 #if defined(WINIFL) || defined(WINPGI)
4403 open (iscpp,file=scpname,status='old',readonly,shared)
4404 #elif (defined CRAY) || (defined AIX)
4405 open (iscpp,file=scpname,status='old',action='read')
4407 open (iscpp,file=scpname,status='old')
4409 open (iscpp,file=scpname,status='old',action='read')
4412 call getenv_loc('PATTERN',patname)
4413 #if defined(WINIFL) || defined(WINPGI)
4414 open (icbase,file=patname,status='old',readonly,shared)
4415 #elif (defined CRAY) || (defined AIX)
4416 open (icbase,file=patname,status='old',action='read')
4418 open (icbase,file=patname,status='old')
4420 open (icbase,file=patname,status='old',action='read')
4423 ! Open output file only for CG processes
4424 ! print *,"Processor",myrank," fg_rank",fg_rank
4425 if (fg_rank.eq.0) then
4427 if (nodes.eq.1) then
4430 npos = dlog10(dfloat(nodes-1))+1
4432 if (npos.lt.3) npos=3
4433 write (liczba,'(i1)') npos
4434 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
4436 write (liczba,form) me
4437 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
4438 liczba(:ilen(liczba))
4439 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4441 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4443 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
4444 liczba(:ilen(liczba))//'.mol2'
4445 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4446 liczba(:ilen(liczba))//'.stat'
4448 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
4449 //liczba(:ilen(liczba))//'.stat')
4450 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
4453 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4454 liczba(:ilen(liczba))//'.const'
4459 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
4460 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
4461 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
4462 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
4463 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
4465 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
4467 rest2name=prefix(:ilen(prefix))//'.rst'
4469 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
4472 #if defined(AIX) || defined(PGI)
4473 if (me.eq.king .or. .not. out1file) &
4474 open(iout,file=outname,status='unknown')
4476 if (fg_rank.gt.0) then
4477 write (liczba,'(i3.3)') myrank/nfgtasks
4478 write (ll,'(bz,i3.3)') fg_rank
4479 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4484 open(igeom,file=intname,status='unknown',position='append')
4485 open(ipdb,file=pdbname,status='unknown')
4486 open(imol2,file=mol2name,status='unknown')
4487 open(istat,file=statname,status='unknown',position='append')
4489 !1out open(iout,file=outname,status='unknown')
4492 if (me.eq.king .or. .not.out1file) &
4493 open(iout,file=outname,status='unknown')
4495 if (fg_rank.gt.0) then
4496 write (liczba,'(i3.3)') myrank/nfgtasks
4497 write (ll,'(bz,i3.3)') fg_rank
4498 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4503 open(igeom,file=intname,status='unknown',access='append')
4504 open(ipdb,file=pdbname,status='unknown')
4505 open(imol2,file=mol2name,status='unknown')
4506 open(istat,file=statname,status='unknown',access='append')
4508 !1out open(iout,file=outname,status='unknown')
4511 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
4512 csa_seed=prefix(:lenpre)//'.CSA.seed'
4513 csa_history=prefix(:lenpre)//'.CSA.history'
4514 csa_bank=prefix(:lenpre)//'.CSA.bank'
4515 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
4516 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
4517 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
4518 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
4519 csa_int=prefix(:lenpre)//'.int'
4520 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
4521 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
4522 csa_in=prefix(:lenpre)//'.CSA.in'
4523 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4526 write (iout,'(80(1h-))')
4527 write (iout,'(30x,a)') "FILE ASSIGNMENT"
4528 write (iout,'(80(1h-))')
4529 write (iout,*) "Input file : ",&
4530 pref_orig(:ilen(pref_orig))//'.inp'
4531 write (iout,*) "Output file : ",&
4532 outname(:ilen(outname))
4534 write (iout,*) "Sidechain potential file : ",&
4535 sidename(:ilen(sidename))
4537 write (iout,*) "SCp potential file : ",&
4538 scpname(:ilen(scpname))
4540 write (iout,*) "Electrostatic potential file : ",&
4541 elename(:ilen(elename))
4542 write (iout,*) "Cumulant coefficient file : ",&
4543 fouriername(:ilen(fouriername))
4544 write (iout,*) "Torsional parameter file : ",&
4545 torname(:ilen(torname))
4546 write (iout,*) "Double torsional parameter file : ",&
4547 tordname(:ilen(tordname))
4548 write (iout,*) "SCCOR parameter file : ",&
4549 sccorname(:ilen(sccorname))
4550 write (iout,*) "Bond & inertia constant file : ",&
4551 bondname(:ilen(bondname))
4552 write (iout,*) "Bending parameter file : ",&
4553 thetname(:ilen(thetname))
4554 write (iout,*) "Rotamer parameter file : ",&
4555 rotname(:ilen(rotname))
4558 write (iout,*) "Thetpdb parameter file : ",&
4559 thetname_pdb(:ilen(thetname_pdb))
4562 write (iout,*) "Threading database : ",&
4563 patname(:ilen(patname))
4565 write (iout,*)" DIRTMP : ",&
4567 write (iout,'(80(1h-))')
4570 end subroutine openunits
4571 !-----------------------------------------------------------------------------
4574 use geometry_data, only: nres,dc
4576 ! implicit real*8 (a-h,o-z)
4577 ! include 'DIMENSIONS'
4578 ! include 'COMMON.CHAIN'
4579 ! include 'COMMON.IOUNITS'
4580 ! include 'COMMON.MD'
4583 ! real(kind=8) :: var,ene
4585 open(irest2,file=rest2name,status='unknown')
4586 read(irest2,*) totT,EK,potE,totE,t_bath
4589 ! AL 4/17/17: Now reading d_t(0,:) too
4591 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4594 ! AL 4/17/17: Now reading d_c(0,:) too
4596 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4599 read (irest2,*) iset
4603 end subroutine readrst
4604 !-----------------------------------------------------------------------------
4605 subroutine read_fragments
4609 use control_data, only:out1file
4612 ! implicit real*8 (a-h,o-z)
4613 ! include 'DIMENSIONS'
4617 ! include 'COMMON.SETUP'
4618 ! include 'COMMON.CHAIN'
4619 ! include 'COMMON.IOUNITS'
4620 ! include 'COMMON.MD'
4621 ! include 'COMMON.CONTROL'
4624 ! real(kind=8) :: var,ene
4626 read(inp,*) nset,nfrag,npair,nfrag_back
4628 !el from module energy
4629 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
4630 if(.not.allocated(wfrag_back)) then
4631 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4632 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4634 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4635 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
4637 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4638 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
4641 if(me.eq.king.or..not.out1file) &
4642 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4643 " nfrag_back",nfrag_back
4645 read(inp,*) mset(iset)
4647 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4649 if(me.eq.king.or..not.out1file) &
4650 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4651 ifrag(2,i,iset), qinfrag(i,iset)
4654 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4656 if(me.eq.king.or..not.out1file) &
4657 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4658 ipair(2,i,iset), qinpair(i,iset)
4661 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4662 wfrag_back(3,i,iset),&
4663 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4664 if(me.eq.king.or..not.out1file) &
4665 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4666 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4670 end subroutine read_fragments
4671 !-----------------------------------------------------------------------------
4673 !-----------------------------------------------------------------------------
4677 ! implicit real*8 (a-h,o-z)
4678 ! include 'DIMENSIONS'
4679 ! include 'COMMON.CSA'
4680 ! include 'COMMON.BANK'
4681 ! include 'COMMON.IOUNITS'
4683 ! integer :: ntf,ik,iw_pdb
4684 ! real(kind=8) :: var,ene
4686 open(icsa_in,file=csa_in,status="old",err=100)
4687 read(icsa_in,*) nconf
4688 read(icsa_in,*) jstart,jend
4689 read(icsa_in,*) nstmax
4690 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4691 read(icsa_in,*) nran0,nran1,irr
4692 read(icsa_in,*) nseed
4693 read(icsa_in,*) ntotal,cut1,cut2
4694 read(icsa_in,*) estop
4695 read(icsa_in,*) icmax,irestart
4696 read(icsa_in,*) ntbankm,dele,difcut
4697 read(icsa_in,*) iref,rmscut,pnccut
4698 read(icsa_in,*) ndiff
4705 end subroutine csa_read
4706 !-----------------------------------------------------------------------------
4707 subroutine initial_write
4710 ! implicit real*8 (a-h,o-z)
4711 ! include 'DIMENSIONS'
4712 ! include 'COMMON.CSA'
4713 ! include 'COMMON.BANK'
4714 ! include 'COMMON.IOUNITS'
4716 ! integer :: ntf,ik,iw_pdb
4717 ! real(kind=8) :: var,ene
4719 open(icsa_seed,file=csa_seed,status="unknown")
4720 write(icsa_seed,*) "seed"
4722 #if defined(AIX) || defined(PGI)
4723 open(icsa_history,file=csa_history,status="unknown",&
4726 open(icsa_history,file=csa_history,status="unknown",&
4729 write(icsa_history,*) nconf
4730 write(icsa_history,*) jstart,jend
4731 write(icsa_history,*) nstmax
4732 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4733 write(icsa_history,*) nran0,nran1,irr
4734 write(icsa_history,*) nseed
4735 write(icsa_history,*) ntotal,cut1,cut2
4736 write(icsa_history,*) estop
4737 write(icsa_history,*) icmax,irestart
4738 write(icsa_history,*) ntbankm,dele,difcut
4739 write(icsa_history,*) iref,rmscut,pnccut
4740 write(icsa_history,*) ndiff
4742 write(icsa_history,*)
4745 open(icsa_bank1,file=csa_bank1,status="unknown")
4746 write(icsa_bank1,*) 0
4750 end subroutine initial_write
4751 !-----------------------------------------------------------------------------
4752 subroutine restart_write
4755 ! implicit real*8 (a-h,o-z)
4756 ! include 'DIMENSIONS'
4757 ! include 'COMMON.IOUNITS'
4758 ! include 'COMMON.CSA'
4759 ! include 'COMMON.BANK'
4761 ! integer :: ntf,ik,iw_pdb
4762 ! real(kind=8) :: var,ene
4764 #if defined(AIX) || defined(PGI)
4765 open(icsa_history,file=csa_history,position="append")
4767 open(icsa_history,file=csa_history,access="append")
4769 write(icsa_history,*)
4770 write(icsa_history,*) "This is restart"
4771 write(icsa_history,*)
4772 write(icsa_history,*) nconf
4773 write(icsa_history,*) jstart,jend
4774 write(icsa_history,*) nstmax
4775 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4776 write(icsa_history,*) nran0,nran1,irr
4777 write(icsa_history,*) nseed
4778 write(icsa_history,*) ntotal,cut1,cut2
4779 write(icsa_history,*) estop
4780 write(icsa_history,*) icmax,irestart
4781 write(icsa_history,*) ntbankm,dele,difcut
4782 write(icsa_history,*) iref,rmscut,pnccut
4783 write(icsa_history,*) ndiff
4784 write(icsa_history,*)
4785 write(icsa_history,*) "irestart is: ", irestart
4787 write(icsa_history,*)
4791 end subroutine restart_write
4792 !-----------------------------------------------------------------------------
4794 !-----------------------------------------------------------------------------
4795 subroutine write_pdb(npdb,titelloc,ee)
4797 ! implicit real*8 (a-h,o-z)
4798 ! include 'DIMENSIONS'
4799 ! include 'COMMON.IOUNITS'
4800 character(len=50) :: titelloc1
4801 character*(*) :: titelloc
4802 character(len=3) :: zahl
4803 character(len=5) :: liczba5
4805 integer :: npdb !,ilen
4809 ! real(kind=8) :: var,ene
4813 if (npdb.lt.1000) then
4814 call numstr(npdb,zahl)
4815 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4817 if (npdb.lt.10000) then
4818 write(liczba5,'(i1,i4)') 0,npdb
4820 write(liczba5,'(i5)') npdb
4822 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4824 call pdbout(ee,titelloc1,ipdb)
4827 end subroutine write_pdb
4828 !-----------------------------------------------------------------------------
4830 !-----------------------------------------------------------------------------
4831 subroutine write_thread_summary
4832 ! Thread the sequence through a database of known structures
4833 use control_data, only: refstr
4835 use energy_data, only: n_ene_comp
4837 ! implicit real*8 (a-h,o-z)
4838 ! include 'DIMENSIONS'
4840 use MPI_data !include 'COMMON.INFO'
4843 ! include 'COMMON.CONTROL'
4844 ! include 'COMMON.CHAIN'
4845 ! include 'COMMON.DBASE'
4846 ! include 'COMMON.INTERACT'
4847 ! include 'COMMON.VAR'
4848 ! include 'COMMON.THREAD'
4849 ! include 'COMMON.FFIELD'
4850 ! include 'COMMON.SBRIDGE'
4851 ! include 'COMMON.HEADER'
4852 ! include 'COMMON.NAMES'
4853 ! include 'COMMON.IOUNITS'
4854 ! include 'COMMON.TIME1'
4856 integer,dimension(maxthread) :: ip
4857 real(kind=8),dimension(0:n_ene) :: energia
4859 integer :: i,j,ii,jj,ipj,ik,kk,ist
4860 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4862 write (iout,'(30x,a/)') &
4863 ' *********** Summary threading statistics ************'
4864 write (iout,'(a)') 'Initial energies:'
4865 write (iout,'(a4,2x,a12,14a14,3a8)') &
4866 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4867 'RMSnat','NatCONT','NNCONT','RMS'
4868 ! Energy sort patterns
4873 enet=ener(n_ene-1,ip(i))
4876 if (ener(n_ene-1,ip(j)).lt.enet) then
4878 enet=ener(n_ene-1,ip(j))
4890 ist=nres_base(2,ii)+ipatt(2,i)
4892 energia(i)=ener0(kk,i)
4894 etot=ener0(n_ene_comp+1,i)
4895 rmsnat=ener0(n_ene_comp+2,i)
4896 rms=ener0(n_ene_comp+3,i)
4897 frac=ener0(n_ene_comp+4,i)
4898 frac_nn=ener0(n_ene_comp+5,i)
4901 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4902 i,str_nam(ii),ist+1,&
4903 (energia(print_order(kk)),kk=1,nprint_ene),&
4904 etot,rmsnat,frac,frac_nn,rms
4906 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
4907 i,str_nam(ii),ist+1,&
4908 (energia(print_order(kk)),kk=1,nprint_ene),etot
4911 write (iout,'(//a)') 'Final energies:'
4912 write (iout,'(a4,2x,a12,17a14,3a8)') &
4913 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
4914 'RMSnat','NatCONT','NNCONT','RMS'
4918 ist=nres_base(2,ii)+ipatt(2,i)
4920 energia(kk)=ener(kk,ik)
4922 etot=ener(n_ene_comp+1,i)
4923 rmsnat=ener(n_ene_comp+2,i)
4924 rms=ener(n_ene_comp+3,i)
4925 frac=ener(n_ene_comp+4,i)
4926 frac_nn=ener(n_ene_comp+5,i)
4927 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4928 i,str_nam(ii),ist+1,&
4929 (energia(print_order(kk)),kk=1,nprint_ene),&
4930 etot,rmsnat,frac,frac_nn,rms
4932 write (iout,'(/a/)') 'IEXAM array:'
4933 write (iout,'(i5)') nexcl
4935 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
4937 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
4938 'Max. time for threading step ',max_time_for_thread,&
4939 'Average time for threading step: ',ave_time_for_thread
4941 end subroutine write_thread_summary
4942 !-----------------------------------------------------------------------------
4943 subroutine write_stat_thread(ithread,ipattern,ist)
4945 use energy_data, only: n_ene_comp
4947 ! implicit real*8 (a-h,o-z)
4948 ! include "DIMENSIONS"
4949 ! include "COMMON.CONTROL"
4950 ! include "COMMON.IOUNITS"
4951 ! include "COMMON.THREAD"
4952 ! include "COMMON.FFIELD"
4953 ! include "COMMON.DBASE"
4954 ! include "COMMON.NAMES"
4955 real(kind=8),dimension(0:n_ene) :: energia
4957 integer :: ithread,ipattern,ist,i
4958 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
4960 #if defined(AIX) || defined(PGI)
4961 open(istat,file=statname,position='append')
4963 open(istat,file=statname,access='append')
4966 energia(i)=ener(i,ithread)
4968 etot=ener(n_ene_comp+1,ithread)
4969 rmsnat=ener(n_ene_comp+2,ithread)
4970 rms=ener(n_ene_comp+3,ithread)
4971 frac=ener(n_ene_comp+4,ithread)
4972 frac_nn=ener(n_ene_comp+5,ithread)
4973 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4974 ithread,str_nam(ipattern),ist+1,&
4975 (energia(print_order(i)),i=1,nprint_ene),&
4976 etot,rmsnat,frac,frac_nn,rms
4979 end subroutine write_stat_thread
4980 !-----------------------------------------------------------------------------
4982 !-----------------------------------------------------------------------------
4983 end module io_config