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 if (counter.eq.5) then
2772 nres_molec(molecule)=nres_molec(molecule)+1
2774 ! sccor(j,iii)=c(j,ires)
2778 ! print *, "ATOM",atom(1:3)
2779 else if (atom.eq."C5'") then
2780 read (card(19:19),'(a1)') sugar
2781 isugar=sugarcode(sugar,ires)
2786 ! print *,ires,istype(ires)
2789 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2792 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2794 ! write (*,*) card(23:27),ires,itype(ires,1)
2795 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2796 atom.ne.'N' .and. atom.ne.'C' .and. &
2797 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2798 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
2799 .and. atom.ne.'P '.and. &
2800 atom(1:1).ne.'H' .and. &
2801 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
2803 ! write (iout,*) "sidechain ",atom
2804 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2805 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
2806 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2808 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2811 else if ((ions).and.(card(1:6).eq.'HETATM')) then
2813 read (card(12:16),*) atom
2814 print *,"HETATOM", atom
2815 read (card(18:20),'(a3)') res
2816 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
2817 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
2818 .or.(atom(1:2).eq.'K ')) &
2821 if (molecule.ne.5) molecprev=molecule
2823 nres_molec(molecule)=nres_molec(molecule)+1
2824 itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
2825 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2829 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
2830 if (ires.eq.0) return
2831 ! Calculate dummy residue coordinates inside the "chain" of a multichain
2834 if (ires_old.ne.ires) nres_molec(molecule)=nres_molec(molecule)-2
2835 ! print *,'I have', nres_molec(:)
2837 do k=1,4 ! ions are without dummy
2838 if (nres_molec(k).eq.0) cycle
2840 ! write (iout,*) i,itype(i,1)
2841 ! if (itype(i,1).eq.ntyp1) then
2842 ! write (iout,*) "dummy",i,itype(i,1)
2844 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
2845 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
2849 if (itype(i,k).eq.ntyp1_molec(k)) then
2850 if (itype(i+1,k).eq.ntyp1_molec(k)) then
2851 if (itype(i+2,k).eq.0) then
2852 ! print *,"masz sieczke"
2854 if (itype(i+2,j).ne.0) then
2856 itype(i+1,j)=ntyp1_molec(j)
2857 nres_molec(k)=nres_molec(k)-1
2858 nres_molec(j)=nres_molec(j)+1
2864 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
2865 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
2866 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
2868 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2869 ! print *,i,'tu dochodze'
2870 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
2878 c(j,i)=c(j,i-1)-1.9d0*e2(j)
2882 dcj=(c(j,i-2)-c(j,i-3))/2.0
2883 if (dcj.eq.0) dcj=1.23591524223
2888 else !itype(i+1,1).eq.ntyp1
2890 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2891 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
2898 c(j,i)=c(j,i+1)-1.9d0*e2(j)
2902 dcj=(c(j,i+3)-c(j,i+2))/2.0
2903 if (dcj.eq.0) dcj=1.23591524223
2908 endif !itype(i+1,1).eq.ntyp1
2909 endif !itype.eq.ntyp1
2913 ! Calculate the CM of the last side chain.
2917 dc(j,ires)=sccor(j,iii)
2920 call sccenter(ires,iii,sccor)
2926 ! print *,"molecule",molecule
2927 if ((itype(nres,1).ne.10)) then
2929 if (molecule.eq.5) molecule=molecprev
2930 itype(nres,molecule)=ntyp1_molec(molecule)
2931 nres_molec(molecule)=nres_molec(molecule)+1
2933 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2934 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
2941 c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
2945 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
2946 c(j,nres)=c(j,nres-1)+dcj
2947 c(j,2*nres)=c(j,nres)
2951 ! print *,'I have',nres, nres_molec(:)
2953 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
2955 if (nres.ne.nres0) then
2956 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
2958 stop "Error nres value in WHAM input"
2961 !---------------------------------
2962 !el reallocate tables
2965 ! hfrag_alloc(j,i)=hfrag(j,i)
2968 ! bfrag_alloc(j,i)=bfrag(j,i)
2974 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
2975 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
2976 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
2977 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
2981 ! hfrag(j,i)=hfrag_alloc(j,i)
2986 ! bfrag(j,i)=bfrag_alloc(j,i)
2989 !el end reallocate tables
2990 !---------------------------------
2998 c(j,2*nres)=c(j,nres)
3001 if (itype(1,1).eq.ntyp1) then
3005 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3006 call refsys(2,3,4,e1,e2,e3,fail)
3013 c(j,1)=c(j,2)-1.9d0*e2(j)
3017 dcj=(c(j,4)-c(j,3))/2.0
3023 ! First lets assign correct dummy to correct type of chain
3025 if (itype(1,1).eq.ntyp1) then
3026 if (itype(2,1).eq.0) then
3028 if (itype(2,j).ne.0) then
3030 itype(1,j)=ntyp1_molec(j)
3031 nres_molec(1)=nres_molec(1)-1
3032 nres_molec(j)=nres_molec(j)+1
3039 print *,'I have',nres, nres_molec(:)
3041 ! Copy the coordinates to reference coordinates
3047 ! Calculate internal coordinates.
3049 write (iout,'(/a)') &
3050 "Cartesian coordinates of the reference structure"
3051 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3052 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3054 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3055 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3056 (c(j,ires+nres),j=1,3)
3059 ! znamy już nres wiec mozna alokowac tablice
3060 ! Calculate internal coordinates.
3061 if(me.eq.king.or..not.out1file)then
3062 write (iout,'(a)') &
3063 "Backbone and SC coordinates as read from the PDB"
3065 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
3066 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3067 (c(j,nres+ires),j=1,3)
3070 ! NOW LETS ROCK! SORTING
3071 allocate(c_temporary(3,2*nres))
3072 allocate(itype_temporary(nres,5))
3073 allocate(molnum(nres))
3074 allocate(istype_temp(nres))
3075 itype_temporary(:,:)=0
3079 if (itype(i,k).ne.0) then
3081 c_temporary(j,seqalingbegin)=c(j,i)
3082 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3085 itype_temporary(seqalingbegin,k)=itype(i,k)
3086 istype_temp(seqalingbegin)=istype(i)
3087 molnum(seqalingbegin)=k
3088 seqalingbegin=seqalingbegin+1
3094 c(j,i)=c_temporary(j,i)
3099 itype(i,k)=itype_temporary(i,k)
3100 istype(i)=istype_temp(i)
3103 if (itype(1,1).eq.ntyp1) then
3107 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3108 call refsys(2,3,4,e1,e2,e3,fail)
3115 c(j,1)=c(j,2)-1.9d0*e2(j)
3119 dcj=(c(j,4)-c(j,3))/2.0
3127 write (iout,'(/a)') &
3128 "Cartesian coordinates of the reference structure after sorting"
3129 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3130 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3132 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3133 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3134 (c(j,ires+nres),j=1,3)
3138 ! print *,seqalingbegin,nres
3139 if(.not.allocated(vbld)) then
3140 allocate(vbld(2*nres))
3145 if(.not.allocated(vbld_inv)) then
3146 allocate(vbld_inv(2*nres))
3152 if(.not.allocated(theta)) then
3153 allocate(theta(nres+2))
3157 if(.not.allocated(phi)) allocate(phi(nres+2))
3158 if(.not.allocated(alph)) allocate(alph(nres+2))
3159 if(.not.allocated(omeg)) allocate(omeg(nres+2))
3160 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3161 if(.not.allocated(phiref)) allocate(phiref(nres+2))
3162 if(.not.allocated(costtab)) allocate(costtab(nres))
3163 if(.not.allocated(sinttab)) allocate(sinttab(nres))
3164 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3165 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3166 if(.not.allocated(xxref)) allocate(xxref(nres))
3167 if(.not.allocated(yyref)) allocate(yyref(nres))
3168 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3169 if(.not.allocated(dc_norm)) then
3170 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3171 allocate(dc_norm(3,0:2*nres+2))
3175 call int_from_cart(.true.,.false.)
3176 call sc_loc_geom(.false.)
3178 thetaref(i)=theta(i)
3188 dc(j,i)=c(j,i+1)-c(j,i)
3189 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3194 dc(j,i+nres)=c(j,i+nres)-c(j,i)
3195 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3197 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3201 ! Copy the coordinates to reference coordinates
3202 ! Splits to single chain if occurs
3206 ! cref(j,i,cou)=c(j,i)
3210 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
3211 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
3212 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3213 !-----------------------------
3217 write (iout,*) "symetr", symetr
3220 !c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3222 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
3225 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3230 cref(j,i,cou)=c(j,i)
3231 cref(j,i+nres,cou)=c(j,i+nres)
3233 chain_rep(j,lll,kkk)=c(j,i)
3234 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
3238 write (iout,*) chain_length
3239 if (chain_length.eq.0) chain_length=nres
3241 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
3242 chain_rep(j,chain_length+nres,symetr) &
3243 =chain_rep(j,chain_length+nres,1)
3246 ! write (iout,*) "spraw lancuchy",chain_length,symetr
3248 ! do kkk=1,chain_length
3249 ! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
3253 ! makes copy of chains
3254 write (iout,*) "symetr", symetr
3259 if (symetr.gt.1) then
3266 write(iout,*) (tabperm(i,kkk),kkk=1,4)
3272 ! write (iout,*) i,icha
3273 do lll=1,chain_length
3275 if (cou.le.nres) then
3277 kupa=mod(lll,chain_length)
3278 iprzes=(kkk-1)*chain_length+lll
3279 if (kupa.eq.0) kupa=chain_length
3280 ! write (iout,*) "kupa", kupa
3281 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
3282 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
3289 !-koniec robienia kopii
3292 write (iout,*) "nowa struktura", nperm
3294 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
3296 cref(3,i,kkk),cref(1,nres+i,kkk),&
3297 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
3299 100 format (//' alpha-carbon coordinates ',&
3300 ' centroid coordinates'/ &
3301 ' ', 6X,'X',11X,'Y',11X,'Z', &
3302 10X,'X',11X,'Y',11X,'Z')
3303 110 format (a,'(',i3,')',6f12.5)
3309 bfrag(i,j)=bfrag(i,j)-ishift
3315 hfrag(i,j)=hfrag(i,j)-ishift
3321 end subroutine readpdb
3322 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3323 !-----------------------------------------------------------------------------
3325 !-----------------------------------------------------------------------------
3326 subroutine read_control
3340 use random, only: random_init
3341 ! implicit real*8 (a-h,o-z)
3342 ! include 'DIMENSIONS'
3344 use prng, only:prng_restart
3346 logical :: OKRandom!, prng_restart
3349 ! include 'COMMON.IOUNITS'
3350 ! include 'COMMON.TIME1'
3351 ! include 'COMMON.THREAD'
3352 ! include 'COMMON.SBRIDGE'
3353 ! include 'COMMON.CONTROL'
3354 ! include 'COMMON.MCM'
3355 ! include 'COMMON.MAP'
3356 ! include 'COMMON.HEADER'
3357 ! include 'COMMON.CSA'
3358 ! include 'COMMON.CHAIN'
3359 ! include 'COMMON.MUCA'
3360 ! include 'COMMON.MD'
3361 ! include 'COMMON.FFIELD'
3362 ! include 'COMMON.INTERACT'
3363 ! include 'COMMON.SETUP'
3364 !el integer :: KDIAG,ICORFL,IXDR
3365 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
3366 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
3367 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
3368 ! character(len=80) :: ucase
3369 character(len=640) :: controlcard
3371 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
3377 read (INP,'(a)') titel
3378 call card_concat(controlcard,.true.)
3379 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
3380 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
3381 call reada(controlcard,'SEED',seed,0.0D0)
3382 call random_init(seed)
3383 ! Set up the time limit (caution! The time must be input in minutes!)
3384 read_cart=index(controlcard,'READ_CART').gt.0
3385 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3386 call readi(controlcard,'SYM',symetr,1)
3387 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
3388 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
3389 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
3390 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
3391 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
3392 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
3393 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
3394 call reada(controlcard,'DRMS',drms,0.1D0)
3395 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3396 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
3397 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
3398 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
3399 write (iout,'(a,f10.1)')'DRMS = ',drms
3400 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
3401 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
3403 call readi(controlcard,'NZ_START',nz_start,0)
3404 call readi(controlcard,'NZ_END',nz_end,0)
3405 call readi(controlcard,'IZ_SC',iz_sc,0)
3406 timlim=60.0D0*timlim
3407 safety = 60.0d0*safety
3410 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3411 !C SHIELD keyword sets if the shielding effect of side-chains is used
3412 !C 0 denots no shielding is used all peptide are equally despite the
3413 !C solvent accesible area
3414 !C 1 the newly introduced function
3415 !C 2 reseved for further possible developement
3416 call readi(controlcard,'SHIELD',shield_mode,0)
3417 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3418 write(iout,*) "shield_mode",shield_mode
3419 !C Varibles set size of box
3420 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
3421 protein=index(controlcard,"PROTEIN").gt.0
3422 ions=index(controlcard,"IONS").gt.0
3423 nucleic=index(controlcard,"NUCLEIC").gt.0
3424 write (iout,*) "with_theta_constr ",with_theta_constr
3425 AFMlog=(index(controlcard,'AFM'))
3426 selfguide=(index(controlcard,'SELFGUIDE'))
3427 print *,'AFMlog',AFMlog,selfguide,"KUPA"
3428 call readi(controlcard,'GENCONSTR',genconstr,0)
3429 call reada(controlcard,'BOXX',boxxsize,100.0d0)
3430 call reada(controlcard,'BOXY',boxysize,100.0d0)
3431 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3432 call readi(controlcard,'TUBEMOD',tubemode,0)
3433 write (iout,*) TUBEmode,"TUBEMODE"
3434 if (TUBEmode.gt.0) then
3435 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3436 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3437 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3438 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3439 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3440 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3441 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3442 buftubebot=bordtubebot+tubebufthick
3443 buftubetop=bordtubetop-tubebufthick
3446 ! CUTOFFF ON ELECTROSTATICS
3447 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
3448 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3449 write(iout,*) "R_CUT_ELE=",r_cut_ele
3450 ! Lipidic parameters
3451 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3452 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3453 if (lipthick.gt.0.0d0) then
3454 bordliptop=(boxzsize+lipthick)/2.0
3455 bordlipbot=bordliptop-lipthick
3456 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3457 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3458 buflipbot=bordlipbot+lipbufthick
3459 bufliptop=bordliptop-lipbufthick
3460 if ((lipbufthick*2.0d0).gt.lipthick) &
3461 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3462 endif !lipthick.gt.0
3463 write(iout,*) "bordliptop=",bordliptop
3464 write(iout,*) "bordlipbot=",bordlipbot
3465 write(iout,*) "bufliptop=",bufliptop
3466 write(iout,*) "buflipbot=",buflipbot
3467 write (iout,*) "SHIELD MODE",shield_mode
3469 !C-------------------------
3470 minim=(index(controlcard,'MINIMIZE').gt.0)
3471 dccart=(index(controlcard,'CART').gt.0)
3472 overlapsc=(index(controlcard,'OVERLAP').gt.0)
3473 overlapsc=.not.overlapsc
3474 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3475 searchsc=.not.searchsc
3476 sideadd=(index(controlcard,'SIDEADD').gt.0)
3477 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3478 outpdb=(index(controlcard,'PDBOUT').gt.0)
3479 outmol2=(index(controlcard,'MOL2OUT').gt.0)
3480 pdbref=(index(controlcard,'PDBREF').gt.0)
3481 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3482 indpdb=index(controlcard,'PDBSTART')
3483 extconf=(index(controlcard,'EXTCONF').gt.0)
3484 call readi(controlcard,'IPRINT',iprint,0)
3485 call readi(controlcard,'MAXGEN',maxgen,10000)
3486 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3487 call readi(controlcard,"KDIAG",kdiag,0)
3488 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3489 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3490 write (iout,*) "RESCALE_MODE",rescale_mode
3491 split_ene=index(controlcard,'SPLIT_ENE').gt.0
3492 if (index(controlcard,'REGULAR').gt.0.0D0) then
3493 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3497 if (index(controlcard,'CHECKGRAD').gt.0) then
3499 if (index(controlcard,'CART').gt.0) then
3501 elseif (index(controlcard,'CARINT').gt.0) then
3506 elseif (index(controlcard,'THREAD').gt.0) then
3508 call readi(controlcard,'THREAD',nthread,0)
3509 if (nthread.gt.0) then
3510 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3513 write (iout,'(a)')'A number has to follow the THREAD keyword.'
3514 stop 'Error termination in Read_Control.'
3516 else if (index(controlcard,'MCMA').gt.0) then
3518 else if (index(controlcard,'MCEE').gt.0) then
3520 else if (index(controlcard,'MULTCONF').gt.0) then
3522 else if (index(controlcard,'MAP').gt.0) then
3524 call readi(controlcard,'MAP',nmap,0)
3525 else if (index(controlcard,'CSA').gt.0) then
3527 !rc else if (index(controlcard,'ZSCORE').gt.0) then
3529 !rc ZSCORE is rm from UNRES, modecalc=9 is available
3532 !fcm else if (index(controlcard,'MCMF').gt.0) then
3534 else if (index(controlcard,'SOFTREG').gt.0) then
3536 else if (index(controlcard,'CHECK_BOND').gt.0) then
3538 else if (index(controlcard,'TEST').gt.0) then
3540 else if (index(controlcard,'MD').gt.0) then
3542 else if (index(controlcard,'RE ').gt.0) then
3546 lmuca=index(controlcard,'MUCA').gt.0
3547 call readi(controlcard,'MUCADYN',mucadyn,0)
3548 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3549 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3551 write (iout,*) 'MUCADYN=',mucadyn
3552 write (iout,*) 'MUCASMOOTH=',muca_smooth
3555 iscode=index(controlcard,'ONE_LETTER')
3556 indphi=index(controlcard,'PHI')
3557 indback=index(controlcard,'BACK')
3558 iranconf=index(controlcard,'RAND_CONF')
3559 i2ndstr=index(controlcard,'USE_SEC_PRED')
3560 gradout=index(controlcard,'GRADOUT').gt.0
3561 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3562 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3563 if (me.eq.king .or. .not.out1file ) &
3564 write (iout,*) "DISTCHAINMAX",distchainmax
3566 if(me.eq.king.or..not.out1file) &
3567 write (iout,'(2a)') diagmeth(kdiag),&
3568 ' routine used to diagonalize matrices.'
3569 if (shield_mode.gt.0) then
3571 !C VSolvSphere the volume of solving sphere
3573 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3574 !C there will be no distinction between proline peptide group and normal peptide
3575 !C group in case of shielding parameters
3576 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3577 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3578 write (iout,*) VSolvSphere,VSolvSphere_div
3579 !C long axis of side chain
3581 long_r_sidechain(i)=vbldsc0(1,i)
3582 short_r_sidechain(i)=sigma0(i)
3583 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3589 end subroutine read_control
3590 !-----------------------------------------------------------------------------
3591 subroutine read_REMDpar
3593 ! Read REMD settings
3600 use control_data, only:out1file
3601 ! implicit real*8 (a-h,o-z)
3602 ! include 'DIMENSIONS'
3603 ! include 'COMMON.IOUNITS'
3604 ! include 'COMMON.TIME1'
3605 ! include 'COMMON.MD'
3608 !el include 'COMMON.LANGEVIN'
3610 !el include 'COMMON.LANGEVIN.lang0'
3612 ! include 'COMMON.INTERACT'
3613 ! include 'COMMON.NAMES'
3614 ! include 'COMMON.GEO'
3615 ! include 'COMMON.REMD'
3616 ! include 'COMMON.CONTROL'
3617 ! include 'COMMON.SETUP'
3618 ! character(len=80) :: ucase
3619 character(len=320) :: controlcard
3620 character(len=3200) :: controlcard1
3621 integer :: iremd_m_total
3624 ! real(kind=8) :: var,ene
3626 if(me.eq.king.or..not.out1file) &
3627 write (iout,*) "REMD setup"
3629 call card_concat(controlcard,.true.)
3630 call readi(controlcard,"NREP",nrep,3)
3631 call readi(controlcard,"NSTEX",nstex,1000)
3632 call reada(controlcard,"RETMIN",retmin,10.0d0)
3633 call reada(controlcard,"RETMAX",retmax,1000.0d0)
3634 mremdsync=(index(controlcard,'SYNC').gt.0)
3635 call readi(controlcard,"NSYN",i_sync_step,100)
3636 restart1file=(index(controlcard,'REST1FILE').gt.0)
3637 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3638 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3639 if(max_cache_traj_use.gt.max_cache_traj) &
3640 max_cache_traj_use=max_cache_traj
3641 if(me.eq.king.or..not.out1file) then
3642 !d if (traj1file) then
3643 !rc caching is in testing - NTWX is not ignored
3644 !d write (iout,*) "NTWX value is ignored"
3645 !d write (iout,*) " trajectory is stored to one file by master"
3646 !d write (iout,*) " before exchange at NSTEX intervals"
3648 write (iout,*) "NREP= ",nrep
3649 write (iout,*) "NSTEX= ",nstex
3650 write (iout,*) "SYNC= ",mremdsync
3651 write (iout,*) "NSYN= ",i_sync_step
3652 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3655 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3656 if (index(controlcard,'TLIST').gt.0) then
3658 call card_concat(controlcard1,.true.)
3659 read(controlcard1,*) (remd_t(i),i=1,nrep)
3660 if(me.eq.king.or..not.out1file) &
3661 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3664 if (index(controlcard,'MLIST').gt.0) then
3666 call card_concat(controlcard1,.true.)
3667 read(controlcard1,*) (remd_m(i),i=1,nrep)
3668 if(me.eq.king.or..not.out1file) then
3669 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3672 iremd_m_total=iremd_m_total+remd_m(i)
3674 write (iout,*) 'Total number of replicas ',iremd_m_total
3677 if(me.eq.king.or..not.out1file) &
3678 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3680 end subroutine read_REMDpar
3681 !-----------------------------------------------------------------------------
3682 subroutine read_MDpar
3686 use control_data, only: r_cut,rlamb,out1file
3688 use geometry_data, only: pi
3690 ! implicit real*8 (a-h,o-z)
3691 ! include 'DIMENSIONS'
3692 ! include 'COMMON.IOUNITS'
3693 ! include 'COMMON.TIME1'
3694 ! include 'COMMON.MD'
3697 !el include 'COMMON.LANGEVIN'
3699 !el include 'COMMON.LANGEVIN.lang0'
3701 ! include 'COMMON.INTERACT'
3702 ! include 'COMMON.NAMES'
3703 ! include 'COMMON.GEO'
3704 ! include 'COMMON.SETUP'
3705 ! include 'COMMON.CONTROL'
3706 ! include 'COMMON.SPLITELE'
3707 ! character(len=80) :: ucase
3708 character(len=320) :: controlcard
3713 call card_concat(controlcard,.true.)
3714 call readi(controlcard,"NSTEP",n_timestep,1000000)
3715 call readi(controlcard,"NTWE",ntwe,100)
3716 call readi(controlcard,"NTWX",ntwx,1000)
3717 call reada(controlcard,"DT",d_time,1.0d-1)
3718 call reada(controlcard,"DVMAX",dvmax,2.0d1)
3719 call reada(controlcard,"DAMAX",damax,1.0d1)
3720 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3721 call readi(controlcard,"LANG",lang,0)
3722 RESPA = index(controlcard,"RESPA") .gt. 0
3723 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3724 ntime_split0=ntime_split
3725 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3726 ntime_split0=ntime_split
3727 call reada(controlcard,"R_CUT",r_cut,2.0d0)
3728 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3729 rest = index(controlcard,"REST").gt.0
3730 tbf = index(controlcard,"TBF").gt.0
3731 usampl = index(controlcard,"USAMPL").gt.0
3732 mdpdb = index(controlcard,"MDPDB").gt.0
3733 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3734 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3735 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3736 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3737 if (count_reset_moment.eq.0) count_reset_moment=1000000000
3738 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3739 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3740 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3741 if (count_reset_vel.eq.0) count_reset_vel=1000000000
3742 large = index(controlcard,"LARGE").gt.0
3743 print_compon = index(controlcard,"PRINT_COMPON").gt.0
3744 rattle = index(controlcard,"RATTLE").gt.0
3745 ! if performing umbrella sampling, fragments constrained are read from the fragment file
3751 if(me.eq.king.or..not.out1file) then
3753 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3755 write (iout,'(a)') "The units are:"
3756 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3757 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3758 " acceleration: angstrom/(48.9 fs)**2"
3759 write (iout,'(a)') "energy: kcal/mol, temperature: K"
3761 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3762 write (iout,'(a60,f10.5,a)') &
3763 "Initial time step of numerical integration:",d_time,&
3765 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3767 write (iout,'(2a,i4,a)') &
3768 "A-MTS algorithm used; initial time step for fast-varying",&
3769 " short-range forces split into",ntime_split," steps."
3770 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3771 r_cut," lambda",rlamb
3773 write (iout,'(2a,f10.5)') &
3774 "Maximum acceleration threshold to reduce the time step",&
3775 "/increase split number:",damax
3776 write (iout,'(2a,f10.5)') &
3777 "Maximum predicted energy drift to reduce the timestep",&
3778 "/increase split number:",edriftmax
3779 write (iout,'(a60,f10.5)') &
3780 "Maximum velocity threshold to reduce velocities:",dvmax
3781 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3782 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3783 if (rattle) write (iout,'(a60)') &
3784 "Rattle algorithm used to constrain the virtual bonds"
3788 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3789 call reada(controlcard,"RWAT",rwat,1.4d0)
3790 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3791 surfarea=index(controlcard,"SURFAREA").gt.0
3792 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3793 if(me.eq.king.or..not.out1file)then
3794 write (iout,'(/a,$)') "Langevin dynamics calculation"
3796 write (iout,'(a/)') &
3797 " with direct integration of Langevin equations"
3798 else if (lang.eq.2) then
3799 write (iout,'(a/)') " with TINKER stochasic MD integrator"
3800 else if (lang.eq.3) then
3801 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3802 else if (lang.eq.4) then
3803 write (iout,'(a/)') " in overdamped mode"
3805 write (iout,'(//a,i5)') &
3806 "=========== ERROR: Unknown Langevin dynamics mode:",lang
3809 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3810 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3811 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3812 write (iout,'(a60,f10.5)') &
3813 "Scaling factor of the friction forces:",scal_fric
3814 if (surfarea) write (iout,'(2a,i10,a)') &
3815 "Friction coefficients will be scaled by solvent-accessible",&
3816 " surface area every",reset_fricmat," steps."
3818 ! Calculate friction coefficients and bounds of stochastic forces
3819 eta=6*pi*cPoise*etawat
3820 if(me.eq.king.or..not.out1file) &
3821 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3823 gamp=scal_fric*(pstok(i)+rwat)*eta
3824 stdfp=dsqrt(2*Rb*t_bath/d_time)
3825 allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
3827 gamsc(i)=scal_fric*(restok(i,1)+rwat)*eta
3828 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
3830 if(me.eq.king.or..not.out1file)then
3831 write (iout,'(/2a/)') &
3832 "Radii of site types and friction coefficients and std's of",&
3833 " stochastic forces of fully exposed sites"
3834 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
3836 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
3837 gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
3841 if(me.eq.king.or..not.out1file)then
3842 write (iout,'(a)') "Berendsen bath calculation"
3843 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3844 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3846 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3847 count_reset_moment," steps"
3849 write (iout,'(a,i10,a)') &
3850 "Velocities will be reset at random every",count_reset_vel,&
3854 if(me.eq.king.or..not.out1file) &
3855 write (iout,'(a31)') "Microcanonical mode calculation"
3857 if(me.eq.king.or..not.out1file)then
3858 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3860 write(iout,*) "MD running with constraints."
3861 write(iout,*) "Equilibration time ", eq_time, " mtus."
3862 write(iout,*) "Constraining ", nfrag," fragments."
3863 write(iout,*) "Length of each fragment, weight and q0:"
3865 write (iout,*) "Set of restraints #",iset
3867 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3868 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3870 write(iout,*) "constraints between ", npair, "fragments."
3871 write(iout,*) "constraint pairs, weights and q0:"
3873 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3874 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3876 write(iout,*) "angle constraints within ", nfrag_back,&
3877 "backbone fragments."
3878 write(iout,*) "fragment, weights:"
3880 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3881 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3882 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3885 iset=mod(kolor,nset)+1
3888 if(me.eq.king.or..not.out1file) &
3889 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3891 end subroutine read_MDpar
3892 !-----------------------------------------------------------------------------
3896 ! implicit real*8 (a-h,o-z)
3897 ! include 'DIMENSIONS'
3898 ! include 'COMMON.MAP'
3899 ! include 'COMMON.IOUNITS'
3900 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3901 character(len=80) :: mapcard !,ucase
3904 ! real(kind=8) :: var,ene
3907 read (inp,'(a)') mapcard
3908 mapcard=ucase(mapcard)
3909 if (index(mapcard,'PHI').gt.0) then
3911 else if (index(mapcard,'THE').gt.0) then
3913 else if (index(mapcard,'ALP').gt.0) then
3915 else if (index(mapcard,'OME').gt.0) then
3918 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
3919 stop 'Error - illegal variable spec in MAP card.'
3921 call readi (mapcard,'RES1',res1(imap),0)
3922 call readi (mapcard,'RES2',res2(imap),0)
3923 if (res1(imap).eq.0) then
3924 res1(imap)=res2(imap)
3925 else if (res2(imap).eq.0) then
3926 res2(imap)=res1(imap)
3928 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
3929 write (iout,'(a)') &
3930 'Error - illegal definition of variable group in MAP.'
3931 stop 'Error - illegal definition of variable group in MAP.'
3933 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
3934 call reada(mapcard,'TO',ang_to(imap),0.0D0)
3935 call readi(mapcard,'NSTEP',nstep(imap),0)
3936 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
3937 write (iout,'(a)') &
3938 'Illegal boundary and/or step size specification in MAP.'
3939 stop 'Illegal boundary and/or step size specification in MAP.'
3943 end subroutine map_read
3944 !-----------------------------------------------------------------------------
3947 use control_data, only: vdisulf
3949 ! implicit real*8 (a-h,o-z)
3950 ! include 'DIMENSIONS'
3951 ! include 'COMMON.IOUNITS'
3952 ! include 'COMMON.GEO'
3953 ! include 'COMMON.CSA'
3954 ! include 'COMMON.BANK'
3955 ! include 'COMMON.CONTROL'
3956 ! character(len=80) :: ucase
3957 character(len=620) :: mcmcard
3959 ! integer :: ntf,ik,iw_pdb
3960 ! real(kind=8) :: var,ene
3962 call card_concat(mcmcard,.true.)
3964 call readi(mcmcard,'NCONF',nconf,50)
3965 call readi(mcmcard,'NADD',nadd,0)
3966 call readi(mcmcard,'JSTART',jstart,1)
3967 call readi(mcmcard,'JEND',jend,1)
3968 call readi(mcmcard,'NSTMAX',nstmax,500000)
3969 call readi(mcmcard,'N0',n0,1)
3970 call readi(mcmcard,'N1',n1,6)
3971 call readi(mcmcard,'N2',n2,4)
3972 call readi(mcmcard,'N3',n3,0)
3973 call readi(mcmcard,'N4',n4,0)
3974 call readi(mcmcard,'N5',n5,0)
3975 call readi(mcmcard,'N6',n6,10)
3976 call readi(mcmcard,'N7',n7,0)
3977 call readi(mcmcard,'N8',n8,0)
3978 call readi(mcmcard,'N9',n9,0)
3979 call readi(mcmcard,'N14',n14,0)
3980 call readi(mcmcard,'N15',n15,0)
3981 call readi(mcmcard,'N16',n16,0)
3982 call readi(mcmcard,'N17',n17,0)
3983 call readi(mcmcard,'N18',n18,0)
3985 vdisulf=(index(mcmcard,'DYNSS').gt.0)
3987 call readi(mcmcard,'NDIFF',ndiff,2)
3988 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
3989 call readi(mcmcard,'IS1',is1,1)
3990 call readi(mcmcard,'IS2',is2,8)
3991 call readi(mcmcard,'NRAN0',nran0,4)
3992 call readi(mcmcard,'NRAN1',nran1,2)
3993 call readi(mcmcard,'IRR',irr,1)
3994 call readi(mcmcard,'NSEED',nseed,20)
3995 call readi(mcmcard,'NTOTAL',ntotal,10000)
3996 call reada(mcmcard,'CUT1',cut1,2.0d0)
3997 call reada(mcmcard,'CUT2',cut2,5.0d0)
3998 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
3999 call readi(mcmcard,'ICMAX',icmax,3)
4000 call readi(mcmcard,'IRESTART',irestart,0)
4001 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
4004 call reada(mcmcard,'DELE',dele,20.0d0)
4005 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4006 call readi(mcmcard,'IREF',iref,0)
4007 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4008 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4009 call readi(mcmcard,'NCONF_IN',nconf_in,0)
4010 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4011 write (iout,*) "NCONF_IN",nconf_in
4013 end subroutine csaread
4014 !-----------------------------------------------------------------------------
4018 use control_data, only: MaxMoveType
4021 ! implicit real*8 (a-h,o-z)
4022 ! include 'DIMENSIONS'
4023 ! include 'COMMON.MCM'
4024 ! include 'COMMON.MCE'
4025 ! include 'COMMON.IOUNITS'
4026 ! character(len=80) :: ucase
4027 character(len=320) :: mcmcard
4030 ! real(kind=8) :: var,ene
4032 call card_concat(mcmcard,.true.)
4033 call readi(mcmcard,'MAXACC',maxacc,100)
4034 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4035 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4036 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4037 call readi(mcmcard,'MAXREPM',maxrepm,200)
4038 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4039 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4040 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4041 call reada(mcmcard,'E_UP',e_up,5.0D0)
4042 call reada(mcmcard,'DELTE',delte,0.1D0)
4043 call readi(mcmcard,'NSWEEP',nsweep,5)
4044 call readi(mcmcard,'NSTEPH',nsteph,0)
4045 call readi(mcmcard,'NSTEPC',nstepc,0)
4046 call reada(mcmcard,'TMIN',tmin,298.0D0)
4047 call reada(mcmcard,'TMAX',tmax,298.0D0)
4048 call readi(mcmcard,'NWINDOW',nwindow,0)
4049 call readi(mcmcard,'PRINT_MC',print_mc,0)
4050 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4051 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4052 ent_read=(index(mcmcard,'ENT_READ').gt.0)
4053 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4054 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4055 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4056 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4057 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4058 if (nwindow.gt.0) then
4059 allocate(winstart(nwindow)) !!el (maxres)
4060 allocate(winend(nwindow)) !!el
4061 allocate(winlen(nwindow)) !!el
4062 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4064 winlen(i)=winend(i)-winstart(i)+1
4067 if (tmax.lt.tmin) tmax=tmin
4068 if (tmax.eq.tmin) then
4072 if (nstepc.gt.0 .and. nsteph.gt.0) then
4073 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4074 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4076 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4077 ! Probabilities of different move types
4078 sumpro_type(0)=0.0D0
4079 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4080 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4081 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4082 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
4083 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4084 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4085 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4087 print *,'i',i,' sumprotype',sumpro_type(i)
4088 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4089 print *,'i',i,' sumprotype',sumpro_type(i)
4092 end subroutine mcmread
4093 !-----------------------------------------------------------------------------
4094 subroutine read_minim
4097 ! implicit real*8 (a-h,o-z)
4098 ! include 'DIMENSIONS'
4099 ! include 'COMMON.MINIM'
4100 ! include 'COMMON.IOUNITS'
4101 ! character(len=80) :: ucase
4102 character(len=320) :: minimcard
4104 ! integer :: ntf,ik,iw_pdb
4105 ! real(kind=8) :: var,ene
4107 call card_concat(minimcard,.true.)
4108 call readi(minimcard,'MAXMIN',maxmin,2000)
4109 call readi(minimcard,'MAXFUN',maxfun,5000)
4110 call readi(minimcard,'MINMIN',minmin,maxmin)
4111 call readi(minimcard,'MINFUN',minfun,maxmin)
4112 call reada(minimcard,'TOLF',tolf,1.0D-2)
4113 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4114 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4115 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4116 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4117 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4118 'Options in energy minimization:'
4119 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4120 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4121 'MinMin:',MinMin,' MinFun:',MinFun,&
4122 ' TolF:',TolF,' RTolF:',RTolF
4124 end subroutine read_minim
4125 !-----------------------------------------------------------------------------
4126 subroutine openunits
4128 use MD_data, only: usampl
4131 use control_data, only:out1file
4132 use control, only: getenv_loc
4133 ! implicit real*8 (a-h,o-z)
4134 ! include 'DIMENSIONS'
4137 character(len=16) :: form,nodename
4138 integer :: nodelen,ierror,npos
4140 ! include 'COMMON.SETUP'
4141 ! include 'COMMON.IOUNITS'
4142 ! include 'COMMON.MD'
4143 ! include 'COMMON.CONTROL'
4144 integer :: lenpre,lenpot,lentmp !,ilen
4146 character(len=3) :: out1file_text !,ucase
4147 character(len=3) :: ll
4150 ! integer :: ntf,ik,iw_pdb
4151 ! real(kind=8) :: var,ene
4153 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4154 call getenv_loc("PREFIX",prefix)
4156 call getenv_loc("POT",pot)
4157 call getenv_loc("DIRTMP",tmpdir)
4158 call getenv_loc("CURDIR",curdir)
4159 call getenv_loc("OUT1FILE",out1file_text)
4160 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4161 out1file_text=ucase(out1file_text)
4162 if (out1file_text(1:1).eq."Y") then
4165 out1file=fg_rank.gt.0
4170 if (lentmp.gt.0) then
4171 write (*,'(80(1h!))')
4172 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
4173 write (*,'(80(1h!))')
4174 write (*,*)"All output files will be on node /tmp directory."
4176 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4177 if (me.eq.king) then
4178 write (*,*) "The master node is ",nodename
4179 else if (fg_rank.eq.0) then
4180 write (*,*) "I am the CG slave node ",nodename
4182 write (*,*) "I am the FG slave node ",nodename
4185 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
4186 lenpre = lentmp+lenpre+1
4188 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
4189 ! Get the names and open the input files
4190 #if defined(WINIFL) || defined(WINPGI)
4191 open(1,file=pref_orig(:ilen(pref_orig))// &
4192 '.inp',status='old',readonly,shared)
4193 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4194 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4195 ! Get parameter filenames and open the parameter files.
4196 call getenv_loc('BONDPAR',bondname)
4197 open (ibond,file=bondname,status='old',readonly,shared)
4198 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4199 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
4200 call getenv_loc('THETPAR',thetname)
4201 open (ithep,file=thetname,status='old',readonly,shared)
4202 call getenv_loc('ROTPAR',rotname)
4203 open (irotam,file=rotname,status='old',readonly,shared)
4204 call getenv_loc('TORPAR',torname)
4205 open (itorp,file=torname,status='old',readonly,shared)
4206 call getenv_loc('TORDPAR',tordname)
4207 open (itordp,file=tordname,status='old',readonly,shared)
4208 call getenv_loc('FOURIER',fouriername)
4209 open (ifourier,file=fouriername,status='old',readonly,shared)
4210 call getenv_loc('ELEPAR',elename)
4211 open (ielep,file=elename,status='old',readonly,shared)
4212 call getenv_loc('SIDEPAR',sidename)
4213 open (isidep,file=sidename,status='old',readonly,shared)
4215 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4216 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
4217 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4218 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
4219 call getenv_loc('TORPAR_NUCL',torname_nucl)
4220 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
4221 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4222 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
4223 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4224 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
4227 #elif (defined CRAY) || (defined AIX)
4228 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4230 ! print *,"Processor",myrank," opened file 1"
4231 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4232 ! print *,"Processor",myrank," opened file 9"
4233 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4234 ! Get parameter filenames and open the parameter files.
4235 call getenv_loc('BONDPAR',bondname)
4236 open (ibond,file=bondname,status='old',action='read')
4237 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4238 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4240 ! print *,"Processor",myrank," opened file IBOND"
4241 call getenv_loc('THETPAR',thetname)
4242 open (ithep,file=thetname,status='old',action='read')
4243 ! print *,"Processor",myrank," opened file ITHEP"
4244 call getenv_loc('ROTPAR',rotname)
4245 open (irotam,file=rotname,status='old',action='read')
4246 ! print *,"Processor",myrank," opened file IROTAM"
4247 call getenv_loc('TORPAR',torname)
4248 open (itorp,file=torname,status='old',action='read')
4249 ! print *,"Processor",myrank," opened file ITORP"
4250 call getenv_loc('TORDPAR',tordname)
4251 open (itordp,file=tordname,status='old',action='read')
4252 ! print *,"Processor",myrank," opened file ITORDP"
4253 call getenv_loc('SCCORPAR',sccorname)
4254 open (isccor,file=sccorname,status='old',action='read')
4255 ! print *,"Processor",myrank," opened file ISCCOR"
4256 call getenv_loc('FOURIER',fouriername)
4257 open (ifourier,file=fouriername,status='old',action='read')
4258 ! print *,"Processor",myrank," opened file IFOURIER"
4259 call getenv_loc('ELEPAR',elename)
4260 open (ielep,file=elename,status='old',action='read')
4261 ! print *,"Processor",myrank," opened file IELEP"
4262 call getenv_loc('SIDEPAR',sidename)
4263 open (isidep,file=sidename,status='old',action='read')
4265 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4266 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4267 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4268 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4269 call getenv_loc('TORPAR_NUCL',torname_nucl)
4270 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4271 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4272 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4273 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4274 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4276 call getenv_loc('LIPTRANPAR',liptranname)
4277 open (iliptranpar,file=liptranname,status='old',action='read')
4278 call getenv_loc('TUBEPAR',tubename)
4279 open (itube,file=tubename,status='old',action='read')
4281 ! print *,"Processor",myrank," opened file ISIDEP"
4282 ! print *,"Processor",myrank," opened parameter files"
4284 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
4285 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4286 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4287 ! Get parameter filenames and open the parameter files.
4288 call getenv_loc('BONDPAR',bondname)
4289 open (ibond,file=bondname,status='old')
4290 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4291 open (ibond_nucl,file=bondname_nucl,status='old')
4293 call getenv_loc('THETPAR',thetname)
4294 open (ithep,file=thetname,status='old')
4295 call getenv_loc('ROTPAR',rotname)
4296 open (irotam,file=rotname,status='old')
4297 call getenv_loc('TORPAR',torname)
4298 open (itorp,file=torname,status='old')
4299 call getenv_loc('TORDPAR',tordname)
4300 open (itordp,file=tordname,status='old')
4301 call getenv_loc('SCCORPAR',sccorname)
4302 open (isccor,file=sccorname,status='old')
4303 call getenv_loc('FOURIER',fouriername)
4304 open (ifourier,file=fouriername,status='old')
4305 call getenv_loc('ELEPAR',elename)
4306 open (ielep,file=elename,status='old')
4307 call getenv_loc('SIDEPAR',sidename)
4308 open (isidep,file=sidename,status='old')
4310 open (ithep_nucl,file=thetname_nucl,status='old')
4311 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4312 open (irotam_nucl,file=rotname_nucl,status='old')
4313 call getenv_loc('TORPAR_NUCL',torname_nucl)
4314 open (itorp_nucl,file=torname_nucl,status='old')
4315 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4316 open (itordp_nucl,file=tordname_nucl,status='old')
4317 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4318 open (isidep_nucl,file=sidename_nucl,status='old')
4320 call getenv_loc('LIPTRANPAR',liptranname)
4321 open (iliptranpar,file=liptranname,status='old')
4322 call getenv_loc('TUBEPAR',tubename)
4323 open (itube,file=tubename,status='old')
4325 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4327 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4328 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4329 ! Get parameter filenames and open the parameter files.
4330 call getenv_loc('BONDPAR',bondname)
4331 open (ibond,file=bondname,status='old',action='read')
4332 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4333 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4334 call getenv_loc('THETPAR',thetname)
4335 open (ithep,file=thetname,status='old',action='read')
4336 call getenv_loc('ROTPAR',rotname)
4337 open (irotam,file=rotname,status='old',action='read')
4338 call getenv_loc('TORPAR',torname)
4339 open (itorp,file=torname,status='old',action='read')
4340 call getenv_loc('TORDPAR',tordname)
4341 open (itordp,file=tordname,status='old',action='read')
4342 call getenv_loc('SCCORPAR',sccorname)
4343 open (isccor,file=sccorname,status='old',action='read')
4345 call getenv_loc('THETPARPDB',thetname_pdb)
4346 print *,"thetname_pdb ",thetname_pdb
4347 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
4348 print *,ithep_pdb," opened"
4350 call getenv_loc('FOURIER',fouriername)
4351 open (ifourier,file=fouriername,status='old',readonly)
4352 call getenv_loc('ELEPAR',elename)
4353 open (ielep,file=elename,status='old',readonly)
4354 call getenv_loc('SIDEPAR',sidename)
4355 open (isidep,file=sidename,status='old',readonly)
4357 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4358 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4359 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4360 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4361 call getenv_loc('TORPAR_NUCL',torname_nucl)
4362 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4363 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4364 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4365 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4366 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4368 call getenv_loc('LIPTRANPAR',liptranname)
4369 open (iliptranpar,file=liptranname,status='old',action='read')
4370 call getenv_loc('TUBEPAR',tubename)
4371 open (itube,file=tubename,status='old',action='read')
4374 call getenv_loc('ROTPARPDB',rotname_pdb)
4375 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
4378 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
4379 #if defined(WINIFL) || defined(WINPGI)
4380 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
4381 #elif (defined CRAY) || (defined AIX)
4382 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4384 open (iscpp_nucl,file=scpname_nucl,status='old')
4386 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4391 ! 8/9/01 In the newest version SCp interaction constants are read from a file
4392 ! Use -DOLDSCP to use hard-coded constants instead.
4394 call getenv_loc('SCPPAR',scpname)
4395 #if defined(WINIFL) || defined(WINPGI)
4396 open (iscpp,file=scpname,status='old',readonly,shared)
4397 #elif (defined CRAY) || (defined AIX)
4398 open (iscpp,file=scpname,status='old',action='read')
4400 open (iscpp,file=scpname,status='old')
4402 open (iscpp,file=scpname,status='old',action='read')
4405 call getenv_loc('PATTERN',patname)
4406 #if defined(WINIFL) || defined(WINPGI)
4407 open (icbase,file=patname,status='old',readonly,shared)
4408 #elif (defined CRAY) || (defined AIX)
4409 open (icbase,file=patname,status='old',action='read')
4411 open (icbase,file=patname,status='old')
4413 open (icbase,file=patname,status='old',action='read')
4416 ! Open output file only for CG processes
4417 ! print *,"Processor",myrank," fg_rank",fg_rank
4418 if (fg_rank.eq.0) then
4420 if (nodes.eq.1) then
4423 npos = dlog10(dfloat(nodes-1))+1
4425 if (npos.lt.3) npos=3
4426 write (liczba,'(i1)') npos
4427 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
4429 write (liczba,form) me
4430 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
4431 liczba(:ilen(liczba))
4432 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4434 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4436 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
4437 liczba(:ilen(liczba))//'.mol2'
4438 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4439 liczba(:ilen(liczba))//'.stat'
4441 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
4442 //liczba(:ilen(liczba))//'.stat')
4443 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
4446 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4447 liczba(:ilen(liczba))//'.const'
4452 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
4453 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
4454 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
4455 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
4456 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
4458 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
4460 rest2name=prefix(:ilen(prefix))//'.rst'
4462 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
4465 #if defined(AIX) || defined(PGI)
4466 if (me.eq.king .or. .not. out1file) &
4467 open(iout,file=outname,status='unknown')
4469 if (fg_rank.gt.0) then
4470 write (liczba,'(i3.3)') myrank/nfgtasks
4471 write (ll,'(bz,i3.3)') fg_rank
4472 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4477 open(igeom,file=intname,status='unknown',position='append')
4478 open(ipdb,file=pdbname,status='unknown')
4479 open(imol2,file=mol2name,status='unknown')
4480 open(istat,file=statname,status='unknown',position='append')
4482 !1out open(iout,file=outname,status='unknown')
4485 if (me.eq.king .or. .not.out1file) &
4486 open(iout,file=outname,status='unknown')
4488 if (fg_rank.gt.0) then
4489 write (liczba,'(i3.3)') myrank/nfgtasks
4490 write (ll,'(bz,i3.3)') fg_rank
4491 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4496 open(igeom,file=intname,status='unknown',access='append')
4497 open(ipdb,file=pdbname,status='unknown')
4498 open(imol2,file=mol2name,status='unknown')
4499 open(istat,file=statname,status='unknown',access='append')
4501 !1out open(iout,file=outname,status='unknown')
4504 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
4505 csa_seed=prefix(:lenpre)//'.CSA.seed'
4506 csa_history=prefix(:lenpre)//'.CSA.history'
4507 csa_bank=prefix(:lenpre)//'.CSA.bank'
4508 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
4509 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
4510 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
4511 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
4512 csa_int=prefix(:lenpre)//'.int'
4513 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
4514 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
4515 csa_in=prefix(:lenpre)//'.CSA.in'
4516 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4519 write (iout,'(80(1h-))')
4520 write (iout,'(30x,a)') "FILE ASSIGNMENT"
4521 write (iout,'(80(1h-))')
4522 write (iout,*) "Input file : ",&
4523 pref_orig(:ilen(pref_orig))//'.inp'
4524 write (iout,*) "Output file : ",&
4525 outname(:ilen(outname))
4527 write (iout,*) "Sidechain potential file : ",&
4528 sidename(:ilen(sidename))
4530 write (iout,*) "SCp potential file : ",&
4531 scpname(:ilen(scpname))
4533 write (iout,*) "Electrostatic potential file : ",&
4534 elename(:ilen(elename))
4535 write (iout,*) "Cumulant coefficient file : ",&
4536 fouriername(:ilen(fouriername))
4537 write (iout,*) "Torsional parameter file : ",&
4538 torname(:ilen(torname))
4539 write (iout,*) "Double torsional parameter file : ",&
4540 tordname(:ilen(tordname))
4541 write (iout,*) "SCCOR parameter file : ",&
4542 sccorname(:ilen(sccorname))
4543 write (iout,*) "Bond & inertia constant file : ",&
4544 bondname(:ilen(bondname))
4545 write (iout,*) "Bending parameter file : ",&
4546 thetname(:ilen(thetname))
4547 write (iout,*) "Rotamer parameter file : ",&
4548 rotname(:ilen(rotname))
4551 write (iout,*) "Thetpdb parameter file : ",&
4552 thetname_pdb(:ilen(thetname_pdb))
4555 write (iout,*) "Threading database : ",&
4556 patname(:ilen(patname))
4558 write (iout,*)" DIRTMP : ",&
4560 write (iout,'(80(1h-))')
4563 end subroutine openunits
4564 !-----------------------------------------------------------------------------
4567 use geometry_data, only: nres,dc
4569 ! implicit real*8 (a-h,o-z)
4570 ! include 'DIMENSIONS'
4571 ! include 'COMMON.CHAIN'
4572 ! include 'COMMON.IOUNITS'
4573 ! include 'COMMON.MD'
4576 ! real(kind=8) :: var,ene
4578 open(irest2,file=rest2name,status='unknown')
4579 read(irest2,*) totT,EK,potE,totE,t_bath
4582 ! AL 4/17/17: Now reading d_t(0,:) too
4584 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4587 ! AL 4/17/17: Now reading d_c(0,:) too
4589 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4592 read (irest2,*) iset
4596 end subroutine readrst
4597 !-----------------------------------------------------------------------------
4598 subroutine read_fragments
4602 use control_data, only:out1file
4605 ! implicit real*8 (a-h,o-z)
4606 ! include 'DIMENSIONS'
4610 ! include 'COMMON.SETUP'
4611 ! include 'COMMON.CHAIN'
4612 ! include 'COMMON.IOUNITS'
4613 ! include 'COMMON.MD'
4614 ! include 'COMMON.CONTROL'
4617 ! real(kind=8) :: var,ene
4619 read(inp,*) nset,nfrag,npair,nfrag_back
4621 !el from module energy
4622 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
4623 if(.not.allocated(wfrag_back)) then
4624 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4625 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4627 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4628 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
4630 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4631 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
4634 if(me.eq.king.or..not.out1file) &
4635 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4636 " nfrag_back",nfrag_back
4638 read(inp,*) mset(iset)
4640 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4642 if(me.eq.king.or..not.out1file) &
4643 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4644 ifrag(2,i,iset), qinfrag(i,iset)
4647 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4649 if(me.eq.king.or..not.out1file) &
4650 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4651 ipair(2,i,iset), qinpair(i,iset)
4654 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4655 wfrag_back(3,i,iset),&
4656 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4657 if(me.eq.king.or..not.out1file) &
4658 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4659 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4663 end subroutine read_fragments
4664 !-----------------------------------------------------------------------------
4666 !-----------------------------------------------------------------------------
4670 ! implicit real*8 (a-h,o-z)
4671 ! include 'DIMENSIONS'
4672 ! include 'COMMON.CSA'
4673 ! include 'COMMON.BANK'
4674 ! include 'COMMON.IOUNITS'
4676 ! integer :: ntf,ik,iw_pdb
4677 ! real(kind=8) :: var,ene
4679 open(icsa_in,file=csa_in,status="old",err=100)
4680 read(icsa_in,*) nconf
4681 read(icsa_in,*) jstart,jend
4682 read(icsa_in,*) nstmax
4683 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4684 read(icsa_in,*) nran0,nran1,irr
4685 read(icsa_in,*) nseed
4686 read(icsa_in,*) ntotal,cut1,cut2
4687 read(icsa_in,*) estop
4688 read(icsa_in,*) icmax,irestart
4689 read(icsa_in,*) ntbankm,dele,difcut
4690 read(icsa_in,*) iref,rmscut,pnccut
4691 read(icsa_in,*) ndiff
4698 end subroutine csa_read
4699 !-----------------------------------------------------------------------------
4700 subroutine initial_write
4703 ! implicit real*8 (a-h,o-z)
4704 ! include 'DIMENSIONS'
4705 ! include 'COMMON.CSA'
4706 ! include 'COMMON.BANK'
4707 ! include 'COMMON.IOUNITS'
4709 ! integer :: ntf,ik,iw_pdb
4710 ! real(kind=8) :: var,ene
4712 open(icsa_seed,file=csa_seed,status="unknown")
4713 write(icsa_seed,*) "seed"
4715 #if defined(AIX) || defined(PGI)
4716 open(icsa_history,file=csa_history,status="unknown",&
4719 open(icsa_history,file=csa_history,status="unknown",&
4722 write(icsa_history,*) nconf
4723 write(icsa_history,*) jstart,jend
4724 write(icsa_history,*) nstmax
4725 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4726 write(icsa_history,*) nran0,nran1,irr
4727 write(icsa_history,*) nseed
4728 write(icsa_history,*) ntotal,cut1,cut2
4729 write(icsa_history,*) estop
4730 write(icsa_history,*) icmax,irestart
4731 write(icsa_history,*) ntbankm,dele,difcut
4732 write(icsa_history,*) iref,rmscut,pnccut
4733 write(icsa_history,*) ndiff
4735 write(icsa_history,*)
4738 open(icsa_bank1,file=csa_bank1,status="unknown")
4739 write(icsa_bank1,*) 0
4743 end subroutine initial_write
4744 !-----------------------------------------------------------------------------
4745 subroutine restart_write
4748 ! implicit real*8 (a-h,o-z)
4749 ! include 'DIMENSIONS'
4750 ! include 'COMMON.IOUNITS'
4751 ! include 'COMMON.CSA'
4752 ! include 'COMMON.BANK'
4754 ! integer :: ntf,ik,iw_pdb
4755 ! real(kind=8) :: var,ene
4757 #if defined(AIX) || defined(PGI)
4758 open(icsa_history,file=csa_history,position="append")
4760 open(icsa_history,file=csa_history,access="append")
4762 write(icsa_history,*)
4763 write(icsa_history,*) "This is restart"
4764 write(icsa_history,*)
4765 write(icsa_history,*) nconf
4766 write(icsa_history,*) jstart,jend
4767 write(icsa_history,*) nstmax
4768 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4769 write(icsa_history,*) nran0,nran1,irr
4770 write(icsa_history,*) nseed
4771 write(icsa_history,*) ntotal,cut1,cut2
4772 write(icsa_history,*) estop
4773 write(icsa_history,*) icmax,irestart
4774 write(icsa_history,*) ntbankm,dele,difcut
4775 write(icsa_history,*) iref,rmscut,pnccut
4776 write(icsa_history,*) ndiff
4777 write(icsa_history,*)
4778 write(icsa_history,*) "irestart is: ", irestart
4780 write(icsa_history,*)
4784 end subroutine restart_write
4785 !-----------------------------------------------------------------------------
4787 !-----------------------------------------------------------------------------
4788 subroutine write_pdb(npdb,titelloc,ee)
4790 ! implicit real*8 (a-h,o-z)
4791 ! include 'DIMENSIONS'
4792 ! include 'COMMON.IOUNITS'
4793 character(len=50) :: titelloc1
4794 character*(*) :: titelloc
4795 character(len=3) :: zahl
4796 character(len=5) :: liczba5
4798 integer :: npdb !,ilen
4802 ! real(kind=8) :: var,ene
4806 if (npdb.lt.1000) then
4807 call numstr(npdb,zahl)
4808 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4810 if (npdb.lt.10000) then
4811 write(liczba5,'(i1,i4)') 0,npdb
4813 write(liczba5,'(i5)') npdb
4815 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4817 call pdbout(ee,titelloc1,ipdb)
4820 end subroutine write_pdb
4821 !-----------------------------------------------------------------------------
4823 !-----------------------------------------------------------------------------
4824 subroutine write_thread_summary
4825 ! Thread the sequence through a database of known structures
4826 use control_data, only: refstr
4828 use energy_data, only: n_ene_comp
4830 ! implicit real*8 (a-h,o-z)
4831 ! include 'DIMENSIONS'
4833 use MPI_data !include 'COMMON.INFO'
4836 ! include 'COMMON.CONTROL'
4837 ! include 'COMMON.CHAIN'
4838 ! include 'COMMON.DBASE'
4839 ! include 'COMMON.INTERACT'
4840 ! include 'COMMON.VAR'
4841 ! include 'COMMON.THREAD'
4842 ! include 'COMMON.FFIELD'
4843 ! include 'COMMON.SBRIDGE'
4844 ! include 'COMMON.HEADER'
4845 ! include 'COMMON.NAMES'
4846 ! include 'COMMON.IOUNITS'
4847 ! include 'COMMON.TIME1'
4849 integer,dimension(maxthread) :: ip
4850 real(kind=8),dimension(0:n_ene) :: energia
4852 integer :: i,j,ii,jj,ipj,ik,kk,ist
4853 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4855 write (iout,'(30x,a/)') &
4856 ' *********** Summary threading statistics ************'
4857 write (iout,'(a)') 'Initial energies:'
4858 write (iout,'(a4,2x,a12,14a14,3a8)') &
4859 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4860 'RMSnat','NatCONT','NNCONT','RMS'
4861 ! Energy sort patterns
4866 enet=ener(n_ene-1,ip(i))
4869 if (ener(n_ene-1,ip(j)).lt.enet) then
4871 enet=ener(n_ene-1,ip(j))
4883 ist=nres_base(2,ii)+ipatt(2,i)
4885 energia(i)=ener0(kk,i)
4887 etot=ener0(n_ene_comp+1,i)
4888 rmsnat=ener0(n_ene_comp+2,i)
4889 rms=ener0(n_ene_comp+3,i)
4890 frac=ener0(n_ene_comp+4,i)
4891 frac_nn=ener0(n_ene_comp+5,i)
4894 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4895 i,str_nam(ii),ist+1,&
4896 (energia(print_order(kk)),kk=1,nprint_ene),&
4897 etot,rmsnat,frac,frac_nn,rms
4899 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
4900 i,str_nam(ii),ist+1,&
4901 (energia(print_order(kk)),kk=1,nprint_ene),etot
4904 write (iout,'(//a)') 'Final energies:'
4905 write (iout,'(a4,2x,a12,17a14,3a8)') &
4906 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
4907 'RMSnat','NatCONT','NNCONT','RMS'
4911 ist=nres_base(2,ii)+ipatt(2,i)
4913 energia(kk)=ener(kk,ik)
4915 etot=ener(n_ene_comp+1,i)
4916 rmsnat=ener(n_ene_comp+2,i)
4917 rms=ener(n_ene_comp+3,i)
4918 frac=ener(n_ene_comp+4,i)
4919 frac_nn=ener(n_ene_comp+5,i)
4920 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4921 i,str_nam(ii),ist+1,&
4922 (energia(print_order(kk)),kk=1,nprint_ene),&
4923 etot,rmsnat,frac,frac_nn,rms
4925 write (iout,'(/a/)') 'IEXAM array:'
4926 write (iout,'(i5)') nexcl
4928 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
4930 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
4931 'Max. time for threading step ',max_time_for_thread,&
4932 'Average time for threading step: ',ave_time_for_thread
4934 end subroutine write_thread_summary
4935 !-----------------------------------------------------------------------------
4936 subroutine write_stat_thread(ithread,ipattern,ist)
4938 use energy_data, only: n_ene_comp
4940 ! implicit real*8 (a-h,o-z)
4941 ! include "DIMENSIONS"
4942 ! include "COMMON.CONTROL"
4943 ! include "COMMON.IOUNITS"
4944 ! include "COMMON.THREAD"
4945 ! include "COMMON.FFIELD"
4946 ! include "COMMON.DBASE"
4947 ! include "COMMON.NAMES"
4948 real(kind=8),dimension(0:n_ene) :: energia
4950 integer :: ithread,ipattern,ist,i
4951 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
4953 #if defined(AIX) || defined(PGI)
4954 open(istat,file=statname,position='append')
4956 open(istat,file=statname,access='append')
4959 energia(i)=ener(i,ithread)
4961 etot=ener(n_ene_comp+1,ithread)
4962 rmsnat=ener(n_ene_comp+2,ithread)
4963 rms=ener(n_ene_comp+3,ithread)
4964 frac=ener(n_ene_comp+4,ithread)
4965 frac_nn=ener(n_ene_comp+5,ithread)
4966 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4967 ithread,str_nam(ipattern),ist+1,&
4968 (energia(print_order(i)),i=1,nprint_ene),&
4969 etot,rmsnat,frac,frac_nn,rms
4972 end subroutine write_stat_thread
4973 !-----------------------------------------------------------------------------
4975 !-----------------------------------------------------------------------------
4976 end module io_config