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(1)
534 write (iout,*) 'FTORS factor =',ftors(1)
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
557 if ( secstruc(i) .eq. 'H') then
558 ! Helix restraints for this residue
561 phi0(ii) = 45.0D0*deg2rad
562 drange(ii)= 5.0D0*deg2rad
563 phibound(1,i) = phi0(ii)-drange(ii)
564 phibound(2,i) = phi0(ii)+drange(ii)
565 else if (secstruc(i) .eq. 'E') then
566 ! strand restraints for this residue
569 phi0(ii) = 180.0D0*deg2rad
570 drange(ii)= 5.0D0*deg2rad
571 phibound(1,i) = phi0(ii)-drange(ii)
572 phibound(2,i) = phi0(ii)+drange(ii)
574 ! no restraints for this residue
575 ndih_nconstr=ndih_nconstr+1
576 idih_nconstr(ndih_nconstr)=i
580 ! deallocate(secstruc)
583 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
584 ! deallocate(secstruc)
587 write(iout,'(A20)')'Error reading FTORS'
588 ! deallocate(secstruc)
590 end subroutine secstrp2dihc
591 !-----------------------------------------------------------------------------
592 subroutine read_secstr_pred(jin,jout,errors)
594 ! implicit real*8 (a-h,o-z)
595 ! INCLUDE 'DIMENSIONS'
596 ! include 'COMMON.IOUNITS'
597 ! include 'COMMON.CHAIN'
598 !el character(len=1),dimension(nres) :: secstruc !(maxres)
599 !el COMMON/SECONDARYS/secstruc
601 character(len=80) :: line,line1 !,ucase
602 logical :: errflag,errors,blankline
605 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
608 read (jin,'(a)') line
609 write (jout,'(2a)') '> ',line(1:78)
611 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
612 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
613 ! end-groups in the input file "*.spred"
616 do while (index(line1,'$END').eq.0)
617 ! Override commented lines.
620 do while (.not.blankline)
622 call mykey(line,line1,ipos,blankline,errflag)
623 if (errflag) write (jout,'(2a)') &
624 'Error when reading sequence in line: ',line
625 errors=errors .or. errflag
626 if (.not. blankline .and. .not. errflag) then
629 !el if (iseq.le.maxres) then
630 if (line1(1:1).eq.'-' ) then
631 secstruc(iseq)=line1(1:1)
632 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
633 ( ucase(line1(1:1)).eq.'H' ) ) then
634 secstruc(iseq)=ucase(line1(1:1))
637 write (jout,1010) line1(1:1), iseq
642 !el write (jout,1000) iseq,maxres
645 do while (ipos1.le.iend)
650 !el if (iseq.le.maxres) then
651 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
652 secstruc(iseq)=line1(ipos1-1:ipos1-1)
653 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
654 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
655 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
658 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
663 !el write (jout,1000) iseq,maxres
670 read (jin,'(a)') line
671 write (jout,'(2a)') '> ',line(1:78)
675 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
677 !d check whether the found length of the chain is correct.
678 length_of_chain=iseq-1
679 if (length_of_chain .ne. nres) then
681 write (jout,'(a,i4,a,i4,a)') &
682 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
683 ,length_of_chain,') does not match with the number of residues (' &
688 1000 format('Error - the number of residues (',i4,&
689 ') has exceeded maximum (',i4,').')
690 1010 format ('Error - unrecognized secondary structure label',a4,&
693 end subroutine read_secstr_pred
695 !-----------------------------------------------------------------------------
697 !-----------------------------------------------------------------------------
702 use control_data, only:maxterm !,maxtor
706 use control, only: getenv_loc
708 ! Read the parameters of the probability distributions of the virtual-bond
709 ! valence angles and the side chains and energy parameters.
711 ! Important! Energy-term weights ARE NOT read here; they are read from the
712 ! main input file instead, because NO defaults have yet been set for these
715 ! implicit real*8 (a-h,o-z)
716 ! include 'DIMENSIONS'
721 ! include 'COMMON.IOUNITS'
722 ! include 'COMMON.CHAIN'
723 ! include 'COMMON.INTERACT'
724 ! include 'COMMON.GEO'
725 ! include 'COMMON.LOCAL'
726 ! include 'COMMON.TORSION'
727 ! include 'COMMON.SCCOR'
728 ! include 'COMMON.SCROT'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.NAMES'
731 ! include 'COMMON.SBRIDGE'
732 ! include 'COMMON.MD'
733 ! include 'COMMON.SETUP'
734 character(len=1) :: t1,t2,t3
735 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
736 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
737 logical :: lprint,LaTeX,SPLIT_FOURIERTOR
738 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
739 real(kind=8),dimension(13) :: buse
740 character(len=3) :: lancuch !,ucase
742 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj
743 integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp
744 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
745 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
746 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
749 integer :: ichir1,ichir2,ijunk
752 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
753 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
756 ! For printing parameters after they are read set the following in the UNRES
759 ! setenv PRINT_PARM YES
761 ! To print parameters in LaTeX format rather than as ASCII tables:
765 call getenv_loc("PRINT_PARM",lancuch)
766 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 call getenv_loc("LATEX",lancuch)
768 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
770 dwa16=2.0d0**(1.0d0/6.0d0)
772 ! Assign virtual-bond length
775 vblinv2=vblinv*vblinv
777 ! Read the virtual-bond parameters, masses, and moments of inertia
778 ! and Stokes' radii of the peptide group and side chains
780 allocate(dsc(ntyp1)) !(ntyp1)
781 allocate(dsc_inv(ntyp1)) !(ntyp1)
782 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
783 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
784 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
785 allocate(nbondterm(ntyp)) !(ntyp)
786 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 allocate(msc(ntyp+1)) !(ntyp+1)
796 allocate(isc(ntyp+1)) !(ntyp+1)
797 allocate(restok(ntyp+1)) !(ntyp+1)
799 read (ibond,*) vbldp0,akp,mp,ip,pstok
802 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
803 dsc(i) = vbldsc0(1,i)
807 dsc_inv(i)=1.0D0/dsc(i)
816 allocate(msc(ntyp+1,5)) !(ntyp+1)
817 allocate(isc(ntyp+1,5)) !(ntyp+1)
818 allocate(restok(ntyp+1,5)) !(ntyp+1)
820 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
822 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
823 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
824 dsc(i) = vbldsc0(1,i)
828 dsc_inv(i)=1.0D0/dsc(i)
832 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
835 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
836 ! dsc(i) = vbldsc0_nucl(1,i)
840 ! dsc_inv(i)=1.0D0/dsc(i)
843 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
844 ! do i=1,ntyp_molec(2)
845 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
846 ! aksc_nucl(j,i),abond0_nucl(j,i),&
847 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
848 ! dsc(i) = vbldsc0(1,i)
852 ! dsc_inv(i)=1.0D0/dsc(i)
857 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
858 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
860 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
862 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
863 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
865 write (iout,'(13x,3f10.5)') &
866 vbldsc0(j,i),aksc(j,i),abond0(j,i)
870 if (.not.allocated(ichargecat)) allocate (ichargecat(ntyp_molec(5)))
871 if (oldion.eq.1) then
873 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
874 print *,msc(i,5),restok(i,5)
878 read (iion,*) ncatprotparm
879 allocate(catprm(ncatprotparm,4))
881 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
885 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
886 !----------------------------------------------------
887 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
888 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
889 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
890 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
891 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
892 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
902 allocate(liptranene(ntyp))
903 !C reading lipid parameters
904 write (iout,*) "iliptranpar",iliptranpar
906 read(iliptranpar,*) pepliptran
909 read(iliptranpar,*) liptranene(i)
910 print *,liptranene(i)
916 ! Read the parameters of the probability distribution/energy expression
917 ! of the virtual-bond valence angles theta
920 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
921 (bthet(j,i,1,1),j=1,2)
922 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
923 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
924 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
928 athet(1,i,1,-1)=athet(1,i,1,1)
929 athet(2,i,1,-1)=athet(2,i,1,1)
930 bthet(1,i,1,-1)=-bthet(1,i,1,1)
931 bthet(2,i,1,-1)=-bthet(2,i,1,1)
932 athet(1,i,-1,1)=-athet(1,i,1,1)
933 athet(2,i,-1,1)=-athet(2,i,1,1)
934 bthet(1,i,-1,1)=bthet(1,i,1,1)
935 bthet(2,i,-1,1)=bthet(2,i,1,1)
939 athet(1,i,-1,-1)=athet(1,-i,1,1)
940 athet(2,i,-1,-1)=-athet(2,-i,1,1)
941 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
942 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
943 athet(1,i,-1,1)=athet(1,-i,1,1)
944 athet(2,i,-1,1)=-athet(2,-i,1,1)
945 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
946 bthet(2,i,-1,1)=bthet(2,-i,1,1)
947 athet(1,i,1,-1)=-athet(1,-i,1,1)
948 athet(2,i,1,-1)=athet(2,-i,1,1)
949 bthet(1,i,1,-1)=bthet(1,-i,1,1)
950 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
955 polthet(j,i)=polthet(j,-i)
958 gthet(j,i)=gthet(j,-i)
966 'Parameters of the virtual-bond valence angles:'
967 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
968 ' ATHETA0 ',' A1 ',' A2 ',&
971 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
972 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
974 write (iout,'(/a/9x,5a/79(1h-))') &
975 'Parameters of the expression for sigma(theta_c):',&
976 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
977 ' ALPH3 ',' SIGMA0C '
979 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
980 (polthet(j,i),j=0,3),sigc0(i)
982 write (iout,'(/a/9x,5a/79(1h-))') &
983 'Parameters of the second gaussian:',&
984 ' THETA0 ',' SIGMA0 ',' G1 ',&
987 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
988 sig0(i),(gthet(j,i),j=1,3)
992 'Parameters of the virtual-bond valence angles:'
993 write (iout,'(/a/9x,5a/79(1h-))') &
994 'Coefficients of expansion',&
995 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
996 ' b1*10^1 ',' b2*10^1 '
998 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
999 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
1000 (10*bthet(j,i,1,1),j=1,2)
1002 write (iout,'(/a/9x,5a/79(1h-))') &
1003 'Parameters of the expression for sigma(theta_c):',&
1004 ' alpha0 ',' alph1 ',' alph2 ',&
1005 ' alhp3 ',' sigma0c '
1007 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1008 (polthet(j,i),j=0,3),sigc0(i)
1010 write (iout,'(/a/9x,5a/79(1h-))') &
1011 'Parameters of the second gaussian:',&
1012 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1015 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1016 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1022 ! Read the parameters of Utheta determined from ab initio surfaces
1023 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1025 IF (tor_mode.eq.0) THEN
1026 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1027 ntheterm3,nsingle,ndouble
1028 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1030 !----------------------------------------------------
1031 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1032 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1033 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1034 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1035 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1036 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1037 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1038 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1039 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1040 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1041 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1042 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1043 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1044 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1045 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1046 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1047 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1048 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1049 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1050 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1051 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1052 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1053 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1054 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1056 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1058 ithetyp(i)=-ithetyp(-i)
1061 aa0thet(:,:,:,:)=0.0d0
1062 aathet(:,:,:,:,:)=0.0d0
1063 bbthet(:,:,:,:,:,:)=0.0d0
1064 ccthet(:,:,:,:,:,:)=0.0d0
1065 ddthet(:,:,:,:,:,:)=0.0d0
1066 eethet(:,:,:,:,:,:)=0.0d0
1067 ffthet(:,:,:,:,:,:,:)=0.0d0
1068 ggthet(:,:,:,:,:,:,:)=0.0d0
1070 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1072 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1073 ! VAR:1=non-glicyne non-proline 2=proline
1074 ! VAR:negative values for D-aminoacid
1076 do j=-nthetyp,nthetyp
1077 do k=-nthetyp,nthetyp
1078 read (ithep,'(6a)',end=111,err=111) res1
1079 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1080 ! VAR: aa0thet is variable describing the average value of Foureir
1081 ! VAR: expansion series
1082 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1083 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1084 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1085 read (ithep,*,end=111,err=111) &
1086 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1087 read (ithep,*,end=111,err=111) &
1088 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1089 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1090 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1091 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1093 read (ithep,*,end=111,err=111) &
1094 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1095 ffthet(lll,llll,ll,i,j,k,iblock),&
1096 ggthet(llll,lll,ll,i,j,k,iblock),&
1097 ggthet(lll,llll,ll,i,j,k,iblock),&
1098 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1103 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1104 ! coefficients of theta-and-gamma-dependent terms are zero.
1105 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1106 ! RECOMENTDED AFTER VERSION 3.3)
1110 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1111 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1113 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1114 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1117 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1119 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1122 ! AND COMMENT THE LOOPS BELOW
1126 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1127 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1129 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1130 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1133 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1135 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1140 ! Substitution for D aminoacids from symmetry.
1143 do j=-nthetyp,nthetyp
1144 do k=-nthetyp,nthetyp
1145 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1147 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1151 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1152 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1153 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1154 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1160 ffthet(llll,lll,ll,i,j,k,iblock)= &
1161 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1162 ffthet(lll,llll,ll,i,j,k,iblock)= &
1163 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1164 ggthet(llll,lll,ll,i,j,k,iblock)= &
1165 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1166 ggthet(lll,llll,ll,i,j,k,iblock)= &
1167 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1176 ! Control printout of the coefficients of virtual-bond-angle potentials
1179 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1184 write (iout,'(//4a)') &
1185 'Type ',onelett(i),onelett(j),onelett(k)
1186 write (iout,'(//a,10x,a)') " l","a[l]"
1187 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1188 write (iout,'(i2,1pe15.5)') &
1189 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1191 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1192 "b",l,"c",l,"d",l,"e",l
1194 write (iout,'(i2,4(1pe15.5))') m,&
1195 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1196 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1200 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1201 "f+",l,"f-",l,"g+",l,"g-",l
1204 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1205 ffthet(n,m,l,i,j,k,iblock),&
1206 ffthet(m,n,l,i,j,k,iblock),&
1207 ggthet(n,m,l,i,j,k,iblock),&
1208 ggthet(m,n,l,i,j,k,iblock)
1219 !C here will be the apropriate recalibrating for D-aminoacid
1220 read (ithep,*,end=121,err=121) nthetyp
1221 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1222 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1223 do i=-nthetyp+1,nthetyp-1
1224 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
1225 do j=0,nbend_kcc_Tb(i)
1226 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
1230 write (iout,'(a)') &
1231 "Parameters of the valence-only potentials"
1232 do i=-nthetyp+1,nthetyp-1
1233 write (iout,'(2a)') "Type ",toronelet(i)
1234 do k=0,nbend_kcc_Tb(i)
1235 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1241 write (2,*) "Start reading THETA_PDB",ithep_pdb
1243 ! write (2,*) 'i=',i
1244 read (ithep_pdb,*,err=111,end=111) &
1245 a0thet(i),(athet(j,i,1,1),j=1,2),&
1246 (bthet(j,i,1,1),j=1,2)
1247 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1248 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1249 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1250 sigc0(i)=sigc0(i)**2
1253 athet(1,i,1,-1)=athet(1,i,1,1)
1254 athet(2,i,1,-1)=athet(2,i,1,1)
1255 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1256 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1257 athet(1,i,-1,1)=-athet(1,i,1,1)
1258 athet(2,i,-1,1)=-athet(2,i,1,1)
1259 bthet(1,i,-1,1)=bthet(1,i,1,1)
1260 bthet(2,i,-1,1)=bthet(2,i,1,1)
1263 a0thet(i)=a0thet(-i)
1264 athet(1,i,-1,-1)=athet(1,-i,1,1)
1265 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1266 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1267 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1268 athet(1,i,-1,1)=athet(1,-i,1,1)
1269 athet(2,i,-1,1)=-athet(2,-i,1,1)
1270 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1271 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1272 athet(1,i,1,-1)=-athet(1,-i,1,1)
1273 athet(2,i,1,-1)=athet(2,-i,1,1)
1274 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1275 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1276 theta0(i)=theta0(-i)
1280 polthet(j,i)=polthet(j,-i)
1283 gthet(j,i)=gthet(j,-i)
1286 write (2,*) "End reading THETA_PDB"
1290 !--------------- Reading theta parameters for nucleic acid-------
1291 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1292 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1293 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1294 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1295 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1296 nthetyp_nucl+1,nthetyp_nucl+1))
1297 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1298 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1299 nthetyp_nucl+1,nthetyp_nucl+1))
1300 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1301 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1302 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1303 nthetyp_nucl+1,nthetyp_nucl+1))
1304 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1305 nthetyp_nucl+1,nthetyp_nucl+1))
1306 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1307 nthetyp_nucl+1,nthetyp_nucl+1))
1308 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1309 nthetyp_nucl+1,nthetyp_nucl+1))
1310 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1311 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1312 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1313 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1314 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1315 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1317 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1318 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1320 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1322 aa0thet_nucl(:,:,:)=0.0d0
1323 aathet_nucl(:,:,:,:)=0.0d0
1324 bbthet_nucl(:,:,:,:,:)=0.0d0
1325 ccthet_nucl(:,:,:,:,:)=0.0d0
1326 ddthet_nucl(:,:,:,:,:)=0.0d0
1327 eethet_nucl(:,:,:,:,:)=0.0d0
1328 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1329 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1334 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1335 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1336 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1337 read (ithep_nucl,*,end=111,err=111) &
1338 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1339 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1340 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1341 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1342 read (ithep_nucl,*,end=111,err=111) &
1343 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1344 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1345 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1350 !-------------------------------------------
1351 allocate(nlob(ntyp1)) !(ntyp1)
1352 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1353 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1354 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1361 gaussc(:,:,:,:)=0.0D0
1365 ! Read the parameters of the probability distribution/energy expression
1366 ! of the side chains.
1369 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1373 dsc_inv(i)=1.0D0/dsc(i)
1384 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1385 ((blower(k,l,1),l=1,k),k=1,3)
1386 censc(1,1,-i)=censc(1,1,i)
1387 censc(2,1,-i)=censc(2,1,i)
1388 censc(3,1,-i)=-censc(3,1,i)
1390 read (irotam,*,end=112,err=112) bsc(j,i)
1391 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1392 ((blower(k,l,j),l=1,k),k=1,3)
1393 censc(1,j,-i)=censc(1,j,i)
1394 censc(2,j,-i)=censc(2,j,i)
1395 censc(3,j,-i)=-censc(3,j,i)
1396 ! BSC is amplitude of Gaussian
1403 akl=akl+blower(k,m,j)*blower(l,m,j)
1407 if (((k.eq.3).and.(l.ne.3)) &
1408 .or.((l.eq.3).and.(k.ne.3))) then
1409 gaussc(k,l,j,-i)=-akl
1410 gaussc(l,k,j,-i)=-akl
1412 gaussc(k,l,j,-i)=akl
1413 gaussc(l,k,j,-i)=akl
1422 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1425 if (nlobi.gt.0) then
1427 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1428 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1429 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1430 'log h',(bsc(j,i),j=1,nlobi)
1431 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1432 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1434 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1435 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1438 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1439 write (iout,'(a,f10.4,4(16x,f10.4))') &
1440 'Center ',(bsc(j,i),j=1,nlobi)
1441 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1450 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1451 ! added by Urszula Kozlowska 07/11/2007
1453 !el Maximum number of SC local term fitting function coefficiants
1454 !el integer,parameter :: maxsccoef=65
1456 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1459 read (irotam,*,end=112,err=112)
1461 read (irotam,*,end=112,err=112)
1464 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1468 !---------reading nucleic acid parameters for rotamers-------------------
1469 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1470 do i=1,ntyp_molec(2)
1471 read (irotam_nucl,*,end=112,err=112)
1473 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1479 write (iout,*) "Base rotamer parameters"
1480 do i=1,ntyp_molec(2)
1481 write (iout,'(a)') restyp(i,2)
1482 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1487 ! Read the parameters of the probability distribution/energy expression
1488 ! of the side chains.
1490 write (2,*) "Start reading ROTAM_PDB"
1492 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1496 dsc_inv(i)=1.0D0/dsc(i)
1507 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1508 ((blower(k,l,1),l=1,k),k=1,3)
1510 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1511 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1512 ((blower(k,l,j),l=1,k),k=1,3)
1519 akl=akl+blower(k,m,j)*blower(l,m,j)
1529 write (2,*) "End reading ROTAM_PDB"
1535 !C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1536 !C interaction energy of the Gly, Ala, and Pro prototypes.
1538 read (ifourier,*) nloctyp
1539 SPLIT_FOURIERTOR = nloctyp.lt.0
1540 nloctyp = iabs(nloctyp)
1541 !C allocate(b1(2,nres)) !(2,-maxtor:maxtor)
1542 !C allocate(b2(2,nres)) !(2,-maxtor:maxtor)
1543 !C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
1544 !C allocate(ctilde(2,2,nres))
1545 !C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
1546 !C allocate(gtb1(2,nres))
1547 !C allocate(gtb2(2,nres))
1548 !C allocate(cc(2,2,nres))
1549 !C allocate(dd(2,2,nres))
1550 !C allocate(ee(2,2,nres))
1551 !C allocate(gtcc(2,2,nres))
1552 !C allocate(gtdd(2,2,nres))
1553 !C allocate(gtee(2,2,nres))
1556 allocate(itype2loc(-ntyp1:ntyp1))
1557 allocate(iloctyp(-nloctyp:nloctyp))
1558 allocate(bnew1(3,2,-nloctyp:nloctyp))
1559 allocate(bnew2(3,2,-nloctyp:nloctyp))
1560 allocate(ccnew(3,2,-nloctyp:nloctyp))
1561 allocate(ddnew(3,2,-nloctyp:nloctyp))
1562 allocate(e0new(3,-nloctyp:nloctyp))
1563 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1564 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1565 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1566 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1567 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1568 allocate(e0newtor(3,-nloctyp:nloctyp))
1569 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1571 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1572 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1573 itype2loc(ntyp1)=nloctyp
1574 iloctyp(nloctyp)=ntyp1
1576 itype2loc(-i)=-itype2loc(i)
1579 allocate(iloctyp(-nloctyp:nloctyp))
1580 allocate(itype2loc(-ntyp1:ntyp1))
1587 iloctyp(-i)=-iloctyp(i)
1589 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1590 !c write (iout,*) "nloctyp",nloctyp,
1591 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1592 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1593 !c write (iout,*) "nloctyp",nloctyp,
1594 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1597 !c write (iout,*) "NEWCORR",i
1598 read (ifourier,*,end=115,err=115)
1601 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
1604 !c write (iout,*) "NEWCORR BNEW1"
1605 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1608 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
1611 !c write (iout,*) "NEWCORR BNEW2"
1612 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1614 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
1615 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
1617 !c write (iout,*) "NEWCORR CCNEW"
1618 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1620 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
1621 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
1623 !c write (iout,*) "NEWCORR DDNEW"
1624 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1628 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
1632 !c write (iout,*) "NEWCORR EENEW1"
1633 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1635 read (ifourier,*,end=115,err=115) e0new(ii,i)
1637 !c write (iout,*) (e0new(ii,i),ii=1,3)
1639 !c write (iout,*) "NEWCORR EENEW"
1642 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1643 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1644 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1645 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1650 bnew1(ii,1,-i)= bnew1(ii,1,i)
1651 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1652 bnew2(ii,1,-i)= bnew2(ii,1,i)
1653 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1656 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1657 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1658 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1659 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1660 ccnew(ii,1,-i)=ccnew(ii,1,i)
1661 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1662 ddnew(ii,1,-i)=ddnew(ii,1,i)
1663 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1665 e0new(1,-i)= e0new(1,i)
1666 e0new(2,-i)=-e0new(2,i)
1667 e0new(3,-i)=-e0new(3,i)
1669 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1670 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1671 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1672 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1676 write (iout,'(a)') "Coefficients of the multibody terms"
1677 do i=-nloctyp+1,nloctyp-1
1678 write (iout,*) "Type: ",onelet(iloctyp(i))
1679 write (iout,*) "Coefficients of the expansion of B1"
1681 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1683 write (iout,*) "Coefficients of the expansion of B2"
1685 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1687 write (iout,*) "Coefficients of the expansion of C"
1688 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1689 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1690 write (iout,*) "Coefficients of the expansion of D"
1691 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1692 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1693 write (iout,*) "Coefficients of the expansion of E"
1694 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1697 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1702 IF (SPLIT_FOURIERTOR) THEN
1704 !c write (iout,*) "NEWCORR TOR",i
1705 read (ifourier,*,end=115,err=115)
1708 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
1711 !c write (iout,*) "NEWCORR BNEW1 TOR"
1712 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1715 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
1718 !c write (iout,*) "NEWCORR BNEW2 TOR"
1719 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1721 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
1722 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
1724 !c write (iout,*) "NEWCORR CCNEW TOR"
1725 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1727 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
1728 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
1730 !c write (iout,*) "NEWCORR DDNEW TOR"
1731 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1735 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
1739 !c write (iout,*) "NEWCORR EENEW1 TOR"
1740 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1742 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
1744 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1746 !c write (iout,*) "NEWCORR EENEW TOR"
1749 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1750 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1751 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1752 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1757 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1758 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1759 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1760 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1763 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1764 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1765 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1766 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1768 e0newtor(1,-i)= e0newtor(1,i)
1769 e0newtor(2,-i)=-e0newtor(2,i)
1770 e0newtor(3,-i)=-e0newtor(3,i)
1772 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1773 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1774 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1775 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1779 write (iout,'(a)') &
1780 "Single-body coefficients of the torsional potentials"
1781 do i=-nloctyp+1,nloctyp-1
1782 write (iout,*) "Type: ",onelet(iloctyp(i))
1783 write (iout,*) "Coefficients of the expansion of B1tor"
1785 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1786 j,(bnew1tor(k,j,i),k=1,3)
1788 write (iout,*) "Coefficients of the expansion of B2tor"
1790 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1791 j,(bnew2tor(k,j,i),k=1,3)
1793 write (iout,*) "Coefficients of the expansion of Ctor"
1794 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1795 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1796 write (iout,*) "Coefficients of the expansion of Dtor"
1797 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1798 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1799 write (iout,*) "Coefficients of the expansion of Etor"
1800 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1803 write (iout,'(1hE,2i1,2f10.5)') &
1804 j,k,(eenewtor(l,j,k,i),l=1,2)
1810 do i=-nloctyp+1,nloctyp-1
1813 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1818 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1822 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1823 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1824 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1825 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1828 ENDIF !SPLIT_FOURIER_TOR
1830 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1831 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1832 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1833 allocate(b(13,-nloctyp-1:nloctyp+1))
1835 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1837 read (ifourier,*,end=115,err=115)
1838 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1840 write (iout,*) 'Type ',onelet(iloctyp(i))
1841 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1855 !c B1tilde(1,i) = b(3)
1856 !c! B1tilde(2,i) =-b(5)
1857 !c! B1tilde(1,-i) =-b(3)
1858 !c! B1tilde(2,-i) =b(5)
1859 !c! b1tilde(1,i)=0.0d0
1860 !c b1tilde(2,i)=0.0d0
1865 !cc B1tilde(1,i) = b(3,i)
1866 !cc B1tilde(2,i) =-b(5,i)
1867 !c B1tilde(1,-i) =-b(3,i)
1868 !c B1tilde(2,-i) =b(5,i)
1869 !cc b1tilde(1,i)=0.0d0
1870 !cc b1tilde(2,i)=0.0d0
1871 !cc B2(1,i) = b(2,i)
1872 !cc B2(2,i) = b(4,i)
1874 !c B2(2,-i) =-b(4,i)
1878 CCold(1,1,i)= b(7,i)
1879 CCold(2,2,i)=-b(7,i)
1880 CCold(2,1,i)= b(9,i)
1881 CCold(1,2,i)= b(9,i)
1882 CCold(1,1,-i)= b(7,i)
1883 CCold(2,2,-i)=-b(7,i)
1884 CCold(2,1,-i)=-b(9,i)
1885 CCold(1,2,-i)=-b(9,i)
1890 !c Ctilde(1,1,i)= CCold(1,1,i)
1891 !c Ctilde(1,2,i)= CCold(1,2,i)
1892 !c Ctilde(2,1,i)=-CCold(2,1,i)
1893 !c Ctilde(2,2,i)=-CCold(2,2,i)
1898 !c Ctilde(1,1,i)= CCold(1,1,i)
1899 !c Ctilde(1,2,i)= CCold(1,2,i)
1900 !c Ctilde(2,1,i)=-CCold(2,1,i)
1901 !c Ctilde(2,2,i)=-CCold(2,2,i)
1903 !c Ctilde(1,1,i)=0.0d0
1904 !c Ctilde(1,2,i)=0.0d0
1905 !c Ctilde(2,1,i)=0.0d0
1906 !c Ctilde(2,2,i)=0.0d0
1907 DDold(1,1,i)= b(6,i)
1908 DDold(2,2,i)=-b(6,i)
1909 DDold(2,1,i)= b(8,i)
1910 DDold(1,2,i)= b(8,i)
1911 DDold(1,1,-i)= b(6,i)
1912 DDold(2,2,-i)=-b(6,i)
1913 DDold(2,1,-i)=-b(8,i)
1914 DDold(1,2,-i)=-b(8,i)
1919 !c Dtilde(1,1,i)= DD(1,1,i)
1920 !c Dtilde(1,2,i)= DD(1,2,i)
1921 !c Dtilde(2,1,i)=-DD(2,1,i)
1922 !c Dtilde(2,2,i)=-DD(2,2,i)
1924 !c Dtilde(1,1,i)=0.0d0
1925 !c Dtilde(1,2,i)=0.0d0
1926 !c Dtilde(2,1,i)=0.0d0
1927 !c Dtilde(2,2,i)=0.0d0
1928 EEold(1,1,i)= b(10,i)+b(11,i)
1929 EEold(2,2,i)=-b(10,i)+b(11,i)
1930 EEold(2,1,i)= b(12,i)-b(13,i)
1931 EEold(1,2,i)= b(12,i)+b(13,i)
1932 EEold(1,1,-i)= b(10,i)+b(11,i)
1933 EEold(2,2,-i)=-b(10,i)+b(11,i)
1934 EEold(2,1,-i)=-b(12,i)+b(13,i)
1935 EEold(1,2,-i)=-b(12,i)-b(13,i)
1936 write(iout,*) "TU DOCHODZE"
1942 !c ee(2,1,i)=ee(1,2,i)
1947 "Coefficients of the cumulants (independent of valence angles)"
1948 do i=-nloctyp+1,nloctyp-1
1949 write (iout,*) 'Type ',onelet(iloctyp(i))
1951 write(iout,'(2f10.5)') B(3,i),B(5,i)
1953 write(iout,'(2f10.5)') B(2,i),B(4,i)
1956 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1960 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1964 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1973 ! Read torsional parameters in old format
1975 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1977 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1978 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1979 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1981 !el from energy module--------
1982 allocate(v1(nterm_old,ntortyp,ntortyp))
1983 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1984 !el---------------------------
1989 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1995 write (iout,'(/a/)') 'Torsional constants:'
1998 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1999 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
2005 ! Read torsional parameters
2007 IF (TOR_MODE.eq.0) THEN
2008 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2009 read (itorp,*,end=113,err=113) ntortyp
2010 !el from energy module---------
2011 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2012 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2014 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2015 allocate(vlor2(maxlor,ntortyp,ntortyp))
2016 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
2017 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2019 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
2020 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2021 !el---------------------------
2024 !el---------------------------
2026 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2028 itortyp(i)=-itortyp(-i)
2030 itortyp(ntyp1)=ntortyp
2031 itortyp(-ntyp1)=-ntortyp
2033 write (iout,*) 'ntortyp',ntortyp
2035 do j=-ntortyp+1,ntortyp-1
2036 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
2038 nterm(-i,-j,iblock)=nterm(i,j,iblock)
2039 nlor(-i,-j,iblock)=nlor(i,j,iblock)
2042 do k=1,nterm(i,j,iblock)
2043 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
2045 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
2046 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
2047 v0ij=v0ij+si*v1(k,i,j,iblock)
2049 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
2050 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
2051 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
2053 do k=1,nlor(i,j,iblock)
2054 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
2055 vlor2(k,i,j),vlor3(k,i,j)
2056 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
2059 v0(-i,-j,iblock)=v0ij
2065 write (iout,'(/a/)') 'Torsional constants:'
2067 do i=-ntortyp,ntortyp
2068 do j=-ntortyp,ntortyp
2069 write (iout,*) 'ityp',i,' jtyp',j
2070 write (iout,*) 'Fourier constants'
2071 do k=1,nterm(i,j,iblock)
2072 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
2075 write (iout,*) 'Lorenz constants'
2076 do k=1,nlor(i,j,iblock)
2077 write (iout,'(3(1pe15.5))') &
2078 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2084 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2086 ! 6/23/01 Read parameters for double torsionals
2088 !el from energy module------------
2089 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2090 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2091 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2092 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2093 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2094 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2095 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2096 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2097 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2098 !---------------------------------
2102 do j=-ntortyp+1,ntortyp-1
2103 do k=-ntortyp+1,ntortyp-1
2104 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
2105 ! write (iout,*) "OK onelett",
2108 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
2109 .or. t3.ne.toronelet(k)) then
2110 write (iout,*) "Error in double torsional parameter file",&
2113 call MPI_Finalize(Ierror)
2115 stop "Error in double torsional parameter file"
2117 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
2118 ntermd_2(i,j,k,iblock)
2119 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
2120 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
2121 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
2122 ntermd_1(i,j,k,iblock))
2123 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
2124 ntermd_1(i,j,k,iblock))
2125 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
2126 ntermd_1(i,j,k,iblock))
2127 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
2128 ntermd_1(i,j,k,iblock))
2129 ! Martix of D parameters for one dimesional foureir series
2130 do l=1,ntermd_1(i,j,k,iblock)
2131 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
2132 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
2133 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
2134 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
2135 ! write(iout,*) "whcodze" ,
2136 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
2138 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
2139 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
2140 v2s(m,l,i,j,k,iblock),&
2141 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
2142 ! Martix of D parameters for two dimesional fourier series
2143 do l=1,ntermd_2(i,j,k,iblock)
2145 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
2146 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
2147 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
2148 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
2157 write (iout,*) 'Constants for double torsionals'
2160 do j=-ntortyp+1,ntortyp-1
2161 do k=-ntortyp+1,ntortyp-1
2162 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
2163 ' nsingle',ntermd_1(i,j,k,iblock),&
2164 ' ndouble',ntermd_2(i,j,k,iblock)
2166 write (iout,*) 'Single angles:'
2167 do l=1,ntermd_1(i,j,k,iblock)
2168 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
2169 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
2170 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
2171 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
2174 write (iout,*) 'Pairs of angles:'
2175 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2176 do l=1,ntermd_2(i,j,k,iblock)
2177 write (iout,'(i5,20f10.5)') &
2178 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
2181 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2182 do l=1,ntermd_2(i,j,k,iblock)
2183 write (iout,'(i5,20f10.5)') &
2184 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
2185 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
2195 itype2loc(i)=itortyp(i)
2199 ELSE IF (TOR_MODE.eq.1) THEN
2201 !C read valence-torsional parameters
2202 read (itorp,*,end=121,err=121) ntortyp
2204 write (iout,*) "Valence-torsional parameters read in ntortyp",&
2206 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
2207 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
2210 itype2loc(i)=itortyp(i)
2214 itortyp(i)=-itortyp(-i)
2216 do i=-ntortyp+1,ntortyp-1
2217 do j=-ntortyp+1,ntortyp-1
2218 !C first we read the cos and sin gamma parameters
2219 read (itorp,'(13x,a)',end=121,err=121) string
2220 write (iout,*) i,j,string
2221 read (itorp,*,end=121,err=121) &
2222 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
2223 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
2224 do k=1,nterm_kcc(j,i)
2225 do l=1,nterm_kcc_Tb(j,i)
2226 do ll=1,nterm_kcc_Tb(j,i)
2227 read (itorp,*,end=121,err=121) ii,jj,kk, &
2228 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2236 !c AL 4/8/16: Calculate coefficients from one-body parameters
2238 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2239 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
2240 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
2241 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2242 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2245 print *,i,itortyp(i)
2246 itortyp(i)=itype2loc(i)
2249 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
2251 do i=-ntortyp+1,ntortyp-1
2252 do j=-ntortyp+1,ntortyp-1
2255 do k=1,nterm_kcc_Tb(j,i)
2256 do l=1,nterm_kcc_Tb(j,i)
2257 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
2258 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
2259 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
2260 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
2263 do k=1,nterm_kcc_Tb(j,i)
2264 do l=1,nterm_kcc_Tb(j,i)
2266 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2267 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2268 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2269 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2271 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2272 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2273 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2274 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2278 !c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
2282 write (iout,*) "TOR_MODE>1 only with NEWCORR"
2287 if (tor_mode.gt.0 .and. lprint) then
2288 !c Print valence-torsional parameters
2289 write (iout,'(a)') &
2290 "Parameters of the valence-torsional potentials"
2291 do i=-ntortyp+1,ntortyp-1
2292 do j=-ntortyp+1,ntortyp-1
2293 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
2294 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
2295 do k=1,nterm_kcc(j,i)
2296 do l=1,nterm_kcc_Tb(j,i)
2297 do ll=1,nterm_kcc_Tb(j,i)
2298 write (iout,'(3i5,2f15.4)')&
2299 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2307 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2308 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
2309 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2310 !el from energy module---------
2311 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2312 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2314 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2315 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2316 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2317 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2319 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2320 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2321 !el---------------------------
2324 !el--------------------
2325 read (itorp_nucl,*,end=113,err=113) &
2326 (itortyp_nucl(i),i=1,ntyp_molec(2))
2327 ! print *,itortyp_nucl(:)
2328 !c write (iout,*) 'ntortyp',ntortyp
2331 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
2332 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2335 do k=1,nterm_nucl(i,j)
2336 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2337 v0ij=v0ij+si*v1_nucl(k,i,j)
2340 do k=1,nlor_nucl(i,j)
2341 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
2342 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2343 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2349 ! Read of Side-chain backbone correlation parameters
2350 ! Modified 11 May 2012 by Adasko
2353 read (isccor,*,end=119,err=119) nsccortyp
2355 !el from module energy-------------
2356 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2357 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2358 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2359 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
2360 !-----------------------------------
2362 !el from module energy-------------
2363 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2365 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2367 isccortyp(i)=-isccortyp(-i)
2369 iscprol=isccortyp(20)
2370 ! write (iout,*) 'ntortyp',ntortyp
2372 !c maxinter is maximum interaction sites
2373 !el from module energy---------
2374 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2375 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2376 -nsccortyp:nsccortyp))
2377 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2378 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2379 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2380 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2381 !-----------------------------------
2385 read (isccor,*,end=119,err=119) &
2386 nterm_sccor(i,j),nlor_sccor(i,j)
2392 nterm_sccor(-i,j)=nterm_sccor(i,j)
2393 nterm_sccor(-i,-j)=nterm_sccor(i,j)
2394 nterm_sccor(i,-j)=nterm_sccor(i,j)
2395 do k=1,nterm_sccor(i,j)
2396 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2398 if (j.eq.iscprol) then
2399 if (i.eq.isccortyp(10)) then
2400 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2401 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2403 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2404 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2405 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2406 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2407 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2408 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2409 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2410 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2413 if (i.eq.isccortyp(10)) then
2414 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2415 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2417 if (j.eq.isccortyp(10)) then
2418 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2419 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2421 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2422 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2423 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2424 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2425 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2426 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2430 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2431 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2432 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2433 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2436 do k=1,nlor_sccor(i,j)
2437 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2438 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2439 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2440 (1+vlor3sccor(k,i,j)**2)
2442 v0sccor(l,i,j)=v0ijsccor
2443 v0sccor(l,-i,j)=v0ijsccor1
2444 v0sccor(l,i,-j)=v0ijsccor2
2445 v0sccor(l,-i,-j)=v0ijsccor3
2451 !el from module energy-------------
2452 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
2454 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2455 ! write (iout,*) 'ntortyp',ntortyp
2457 !c maxinter is maximum interaction sites
2458 !el from module energy---------
2459 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2460 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
2461 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2462 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2463 !-----------------------------------
2467 read (isccor,*,end=119,err=119) &
2468 nterm_sccor(i,j),nlor_sccor(i,j)
2472 do k=1,nterm_sccor(i,j)
2473 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2475 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2478 do k=1,nlor_sccor(i,j)
2479 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2480 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2481 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2482 (1+vlor3sccor(k,i,j)**2)
2484 v0sccor(l,i,j)=v0ijsccor !el ,iblock
2493 write (iout,'(/a/)') 'Torsional constants:'
2496 write (iout,*) 'ityp',i,' jtyp',j
2497 write (iout,*) 'Fourier constants'
2498 do k=1,nterm_sccor(i,j)
2499 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
2501 write (iout,*) 'Lorenz constants'
2502 do k=1,nlor_sccor(i,j)
2503 write (iout,'(3(1pe15.5))') &
2504 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2511 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
2512 ! interaction energy of the Gly, Ala, and Pro prototypes.
2515 ! Read electrostatic-interaction parameters
2520 write (iout,'(/a)') 'Electrostatic interaction constants:'
2521 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2522 'IT','JT','APP','BPP','AEL6','AEL3'
2524 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2525 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2526 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2527 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2532 app (i,j)=epp(i,j)*rri*rri
2533 bpp (i,j)=-2.0D0*epp(i,j)*rri
2534 ael6(i,j)=elpp6(i,j)*4.2D0**6
2535 ael3(i,j)=elpp3(i,j)*4.2D0**3
2537 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2543 ! Read side-chain interaction parameters.
2545 !el from module energy - COMMON.INTERACT-------
2546 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2547 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2548 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2550 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2551 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2552 allocate(epslip(ntyp,ntyp))
2560 !--------------------------------
2562 read (isidep,*,end=117,err=117) ipot,expon
2563 if (ipot.lt.1 .or. ipot.gt.5) then
2564 write (iout,'(2a)') 'Error while reading SC interaction',&
2565 'potential file - unknown potential type.'
2567 call MPI_Finalize(Ierror)
2573 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2574 ', exponents are ',expon,2*expon
2575 ! goto (10,20,30,30,40) ipot
2577 !----------------------- LJ potential ---------------------------------
2579 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2580 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2581 (sigma0(i),i=1,ntyp)
2583 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2584 write (iout,'(a/)') 'The epsilon array:'
2585 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2586 write (iout,'(/a)') 'One-body parameters:'
2587 write (iout,'(a,4x,a)') 'residue','sigma'
2588 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2591 !----------------------- LJK potential --------------------------------
2593 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2594 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2595 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2597 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2598 write (iout,'(a/)') 'The epsilon array:'
2599 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2600 write (iout,'(/a)') 'One-body parameters:'
2601 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2602 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2606 !---------------------- GB or BP potential -----------------------------
2609 ! print *,"I AM in SCELE",scelemode
2610 if (scelemode.eq.0) then
2612 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2614 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2615 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2616 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2617 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2619 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2622 ! For the GB potential convert sigma'**2 into chi'
2625 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2629 write (iout,'(/a/)') 'Parameters of the BP potential:'
2630 write (iout,'(a/)') 'The epsilon array:'
2631 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2632 write (iout,'(/a)') 'One-body parameters:'
2633 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2635 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2636 chip(i),alp(i),i=1,ntyp)
2639 ! print *,ntyp,"NTYP"
2640 allocate(icharge(ntyp1))
2641 ! print *,ntyp,icharge(i)
2643 read (isidep,*) (icharge(i),i=1,ntyp)
2644 print *,ntyp,icharge(i)
2645 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2646 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2647 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2648 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2649 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2650 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2651 allocate(epsintab(ntyp,ntyp))
2652 allocate(dtail(2,ntyp,ntyp))
2653 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2654 allocate(wqdip(2,ntyp,ntyp))
2655 allocate(wstate(4,ntyp,ntyp))
2656 allocate(dhead(2,2,ntyp,ntyp))
2657 allocate(nstate(ntyp,ntyp))
2658 allocate(debaykap(ntyp,ntyp))
2660 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2661 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2665 ! write (*,*) "Im in ALAB", i, " ", j
2667 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2668 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2669 chis(i,j),chis(j,i), & !2 w tej linii
2670 nstate(i,j),(wstate(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
2671 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& ! 4 w tej linii
2672 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2673 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2674 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & ! 5 w tej linii
2675 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2676 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2677 IF ((LaTeX).and.(i.gt.24)) then
2678 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2679 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2680 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2681 chis(i,j),chis(j,i) !2 w tej linii
2683 print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i), wqdip(1,i,j)
2688 IF ((LaTeX).and.(i.gt.24)) then
2689 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2690 dhead(1,1,i,j),dhead(2,1,i,j),& ! 2 w tej linii
2691 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2692 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2693 rborn(i,j),rborn(j,i), & ! 3 w tej linii
2694 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2695 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2702 sigma(i,j) = sigma(j,i)
2703 sigmap1(i,j)=sigmap1(j,i)
2704 sigmap2(i,j)=sigmap2(j,i)
2705 sigiso1(i,j)=sigiso1(j,i)
2706 sigiso2(i,j)=sigiso2(j,i)
2707 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2708 nstate(i,j) = nstate(j,i)
2709 dtail(1,i,j) = dtail(1,j,i)
2710 dtail(2,i,j) = dtail(2,j,i)
2712 alphasur(k,i,j) = alphasur(k,j,i)
2713 wstate(k,i,j) = wstate(k,j,i)
2714 alphiso(k,i,j) = alphiso(k,j,i)
2717 dhead(2,1,i,j) = dhead(1,1,j,i)
2718 dhead(2,2,i,j) = dhead(1,2,j,i)
2719 dhead(1,1,i,j) = dhead(2,1,j,i)
2720 dhead(1,2,i,j) = dhead(2,2,j,i)
2722 epshead(i,j) = epshead(j,i)
2723 sig0head(i,j) = sig0head(j,i)
2726 wqdip(k,i,j) = wqdip(k,j,i)
2729 wquad(i,j) = wquad(j,i)
2730 epsintab(i,j) = epsintab(j,i)
2731 debaykap(i,j)=debaykap(j,i)
2732 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2739 !--------------------- GBV potential -----------------------------------
2741 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2742 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2743 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2744 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2746 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2747 write (iout,'(a/)') 'The epsilon array:'
2748 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2749 write (iout,'(/a)') 'One-body parameters:'
2750 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2751 's||/s_|_^2',' chip ',' alph '
2752 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2753 sigii(i),chip(i),alp(i),i=1,ntyp)
2756 write(iout,*)"Wrong ipot"
2762 !-----------------------------------------------------------------------
2763 ! Calculate the "working" parameters of SC interactions.
2765 !el from module energy - COMMON.INTERACT-------
2766 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2767 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2768 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2769 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2770 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2771 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2773 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2779 if (scelemode.eq.0) then
2788 sc_aa_tube_par(:)=0.0d0
2789 sc_bb_tube_par(:)=0.0d0
2791 !--------------------------------
2796 epslip(i,j)=epslip(j,i)
2799 if (scelemode.eq.0) then
2802 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2803 sigma(j,i)=sigma(i,j)
2804 rs0(i,j)=dwa16*sigma(i,j)
2809 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2810 'Working parameters of the SC interactions:',&
2811 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2816 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2818 ! print *,"SIGMA ZLA?",sigma(i,j)
2826 sigeps=dsign(1.0D0,epsij)
2828 aa_aq(i,j)=epsij*rrij*rrij
2829 ! print *,"ADASKO",epsij,rrij,expon
2830 bb_aq(i,j)=-sigeps*epsij*rrij
2831 aa_aq(j,i)=aa_aq(i,j)
2832 bb_aq(j,i)=bb_aq(i,j)
2833 epsijlip=epslip(i,j)
2834 sigeps=dsign(1.0D0,epsijlip)
2835 epsijlip=dabs(epsijlip)
2836 aa_lip(i,j)=epsijlip*rrij*rrij
2837 bb_lip(i,j)=-sigeps*epsijlip*rrij
2838 aa_lip(j,i)=aa_lip(i,j)
2839 bb_lip(j,i)=bb_lip(i,j)
2840 !C write(iout,*) aa_lip
2841 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2842 sigt1sq=sigma0(i)**2
2843 sigt2sq=sigma0(j)**2
2846 ratsig1=sigt2sq/sigt1sq
2847 ratsig2=1.0D0/ratsig1
2848 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2849 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2850 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2854 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2855 sigmaii(i,j)=rsum_max
2856 sigmaii(j,i)=rsum_max
2858 ! sigmaii(i,j)=r0(i,j)
2859 ! sigmaii(j,i)=r0(i,j)
2861 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2862 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2863 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2864 augm(i,j)=epsij*r_augm**(2*expon)
2865 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2872 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2873 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2874 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2879 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2880 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2881 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2882 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2883 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2884 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2885 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2886 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2887 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2888 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2889 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2890 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2891 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2892 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2893 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2894 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2903 read (isidep_nucl,*) ipot_nucl
2904 ! print *,"TU?!",ipot_nucl
2905 if (ipot_nucl.eq.1) then
2906 do i=1,ntyp_molec(2)
2907 do j=i,ntyp_molec(2)
2908 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2909 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2913 do i=1,ntyp_molec(2)
2914 do j=i,ntyp_molec(2)
2915 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2916 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2917 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2921 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2922 do i=1,ntyp_molec(2)
2923 do j=i,ntyp_molec(2)
2924 rrij=sigma_nucl(i,j)
2928 epsij=4*eps_nucl(i,j)
2929 sigeps=dsign(1.0D0,epsij)
2931 aa_nucl(i,j)=epsij*rrij*rrij
2932 bb_nucl(i,j)=-sigeps*epsij*rrij
2933 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2934 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2935 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2936 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2937 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2938 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2941 aa_nucl(i,j)=aa_nucl(j,i)
2942 bb_nucl(i,j)=bb_nucl(j,i)
2943 ael3_nucl(i,j)=ael3_nucl(j,i)
2944 ael6_nucl(i,j)=ael6_nucl(j,i)
2945 ael63_nucl(i,j)=ael63_nucl(j,i)
2946 ael32_nucl(i,j)=ael32_nucl(j,i)
2947 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2948 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2949 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2950 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2951 eps_nucl(i,j)=eps_nucl(j,i)
2952 sigma_nucl(i,j)=sigma_nucl(j,i)
2953 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2957 write(iout,*) "tube param"
2958 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2959 ccavtubpep,dcavtubpep,tubetranenepep
2960 sigmapeptube=sigmapeptube**6
2961 sigeps=dsign(1.0D0,epspeptube)
2962 epspeptube=dabs(epspeptube)
2963 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2964 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2965 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2967 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2968 ccavtub(i),dcavtub(i),tubetranene(i)
2969 sigmasctube=sigmasctube**6
2970 sigeps=dsign(1.0D0,epssctube)
2971 epssctube=dabs(epssctube)
2972 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2973 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2974 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2976 !-----------------READING SC BASE POTENTIALS-----------------------------
2977 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2978 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2979 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2980 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2981 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2982 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2983 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2984 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2985 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2986 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2987 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2988 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2989 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2990 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2991 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2992 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2995 do i=1,ntyp_molec(1)
2996 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2997 write (*,*) "Im in ", i, " ", j
2998 read(isidep_scbase,*) &
2999 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3000 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3001 write(*,*) "eps",eps_scbase(i,j)
3002 read(isidep_scbase,*) &
3003 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3004 chis_scbase(i,j,1),chis_scbase(i,j,2)
3005 read(isidep_scbase,*) &
3006 dhead_scbasei(i,j), &
3007 dhead_scbasej(i,j), &
3008 rborn_scbasei(i,j),rborn_scbasej(i,j)
3009 read(isidep_scbase,*) &
3010 (wdipdip_scbase(k,i,j),k=1,3), &
3011 (wqdip_scbase(k,i,j),k=1,2)
3012 read(isidep_scbase,*) &
3013 alphapol_scbase(i,j), &
3014 epsintab_scbase(i,j)
3017 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
3018 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
3020 do i=1,ntyp_molec(1)
3021 do j=1,ntyp_molec(2)-1
3022 epsij=eps_scbase(i,j)
3023 rrij=sigma_scbase(i,j)
3028 sigeps=dsign(1.0D0,epsij)
3030 aa_scbase(i,j)=epsij*rrij*rrij
3031 bb_scbase(i,j)=-sigeps*epsij*rrij
3034 !-----------------READING PEP BASE POTENTIALS-------------------
3035 allocate(eps_pepbase(ntyp_molec(2)))
3036 allocate(sigma_pepbase(ntyp_molec(2)))
3037 allocate(chi_pepbase(ntyp_molec(2),2))
3038 allocate(chipp_pepbase(ntyp_molec(2),2))
3039 allocate(alphasur_pepbase(4,ntyp_molec(2)))
3040 allocate(sigmap1_pepbase(ntyp_molec(2)))
3041 allocate(sigmap2_pepbase(ntyp_molec(2)))
3042 allocate(chis_pepbase(ntyp_molec(2),2))
3043 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
3046 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
3047 write (*,*) "Im in ", i, " ", j
3048 read(isidep_pepbase,*) &
3049 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3050 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3051 write(*,*) "eps",eps_pepbase(j)
3052 read(isidep_pepbase,*) &
3053 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3054 chis_pepbase(j,1),chis_pepbase(j,2)
3055 read(isidep_pepbase,*) &
3056 (wdipdip_pepbase(k,j),k=1,3)
3058 allocate(aa_pepbase(ntyp_molec(2)))
3059 allocate(bb_pepbase(ntyp_molec(2)))
3061 do j=1,ntyp_molec(2)-1
3062 epsij=eps_pepbase(j)
3063 rrij=sigma_pepbase(j)
3068 sigeps=dsign(1.0D0,epsij)
3070 aa_pepbase(j)=epsij*rrij*rrij
3071 bb_pepbase(j)=-sigeps*epsij*rrij
3073 !--------------READING SC PHOSPHATE-------------------------------------
3074 allocate(eps_scpho(ntyp_molec(1)))
3075 allocate(sigma_scpho(ntyp_molec(1)))
3076 allocate(chi_scpho(ntyp_molec(1),2))
3077 allocate(chipp_scpho(ntyp_molec(1),2))
3078 allocate(alphasur_scpho(4,ntyp_molec(1)))
3079 allocate(sigmap1_scpho(ntyp_molec(1)))
3080 allocate(sigmap2_scpho(ntyp_molec(1)))
3081 allocate(chis_scpho(ntyp_molec(1),2))
3082 allocate(wqq_scpho(ntyp_molec(1)))
3083 allocate(wqdip_scpho(2,ntyp_molec(1)))
3084 allocate(alphapol_scpho(ntyp_molec(1)))
3085 allocate(epsintab_scpho(ntyp_molec(1)))
3086 allocate(dhead_scphoi(ntyp_molec(1)))
3087 allocate(rborn_scphoi(ntyp_molec(1)))
3088 allocate(rborn_scphoj(ntyp_molec(1)))
3089 allocate(alphi_scpho(ntyp_molec(1)))
3093 do j=1,ntyp_molec(1) ! without U then we will take T for U
3094 write (*,*) "Im in scpho ", i, " ", j
3095 read(isidep_scpho,*) &
3096 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
3097 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
3098 write(*,*) "eps",eps_scpho(j)
3099 read(isidep_scpho,*) &
3100 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
3101 chis_scpho(j,1),chis_scpho(j,2)
3102 read(isidep_scpho,*) &
3103 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
3104 read(isidep_scpho,*) &
3105 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
3109 allocate(aa_scpho(ntyp_molec(1)))
3110 allocate(bb_scpho(ntyp_molec(1)))
3112 do j=1,ntyp_molec(1)
3119 sigeps=dsign(1.0D0,epsij)
3121 aa_scpho(j)=epsij*rrij*rrij
3122 bb_scpho(j)=-sigeps*epsij*rrij
3126 read(isidep_peppho,*) &
3127 eps_peppho,sigma_peppho
3128 read(isidep_peppho,*) &
3129 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
3130 read(isidep_peppho,*) &
3131 (wqdip_peppho(k),k=1,2)
3139 sigeps=dsign(1.0D0,epsij)
3141 aa_peppho=epsij*rrij*rrij
3142 bb_peppho=-sigeps*epsij*rrij
3145 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
3150 ! Define the SC-p interaction constants (hard-coded; old style)
3153 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
3155 ! aad(i,1)=0.3D0*4.0D0**12
3156 ! Following line for constants currently implemented
3157 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
3158 aad(i,1)=1.5D0*4.0D0**12
3159 ! aad(i,1)=0.17D0*5.6D0**12
3161 ! "Soft" SC-p repulsion
3163 ! Following line for constants currently implemented
3164 ! aad(i,1)=0.3D0*4.0D0**6
3165 ! "Hard" SC-p repulsion
3166 bad(i,1)=3.0D0*4.0D0**6
3167 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
3176 ! 8/9/01 Read the SC-p interaction constants from file
3179 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
3182 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
3183 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
3184 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
3185 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
3189 write (iout,*) "Parameters of SC-p interactions:"
3191 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
3192 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
3197 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
3199 do i=1,ntyp_molec(2)
3200 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
3202 do i=1,ntyp_molec(2)
3203 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
3204 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
3206 r0pp=1.12246204830937298142*5.16158
3212 ! Define the constants of the disulfide bridge
3216 ! Old arbitrary potential - commented out.
3221 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
3222 ! energy surface of diethyl disulfide.
3223 ! A. Liwo and U. Kozlowska, 11/24/03
3241 allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
3242 allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
3243 allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
3244 allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
3245 allocate(epsintabcat(ntyp,ntyp))
3246 allocate(dtailcat(2,ntyp,ntyp))
3247 allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
3248 allocate(wqdipcat(2,ntyp,ntyp))
3249 allocate(wstatecat(4,ntyp,ntyp))
3250 allocate(dheadcat(2,2,ntyp,ntyp))
3251 allocate(nstatecat(ntyp,ntyp))
3252 allocate(debaykapcat(ntyp,ntyp))
3254 if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
3255 if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
3256 ! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
3257 if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3258 if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3261 if (.not.allocated(ichargecat)) allocate (ichargecat(ntyp_molec(5)))
3262 ! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
3263 if (oldion.eq.0) then
3264 if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
3265 allocate(icharge(1:ntyp1))
3266 read(iion,*) (icharge(i),i=1,ntyp)
3271 do i=1,ntyp_molec(5)
3272 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
3273 print *,msc(i,5),restok(i,5)
3277 do j=1,ntyp_molec(5)
3279 ! do j=1,ntyp_molec(5)
3280 ! write (*,*) "Im in ALAB", i, " ", j
3282 epscat(i,j),sigmacat(i,j), &
3283 ! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
3284 chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
3286 (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
3287 ! chiscat(i,j),chiscat(j,i), &
3288 chis1cat(i,j),chis2cat(i,j), &
3290 nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
3291 dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
3292 dtailcat(1,i,j),dtailcat(2,i,j), &
3293 epsheadcat(i,j),sig0headcat(i,j), &
3295 ! rborncat(i,j),rborncat(j,i),&
3296 rborn1cat(i,j),rborn2cat(i,j),&
3297 (wqdipcat(k,i,j),k=1,2), &
3298 alphapolcat(i,j),alphapolcat(j,i), &
3299 (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
3300 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
3302 ! write (iout,*) 'i= ', i, ' j= ', j
3303 ! write (iout,*) 'epsi0= ', epscat(1,j)
3304 ! write (iout,*) 'sigma0= ', sigmacat(1,j)
3305 ! write (iout,*) 'chi(i,j)= ', chicat(1,j)
3306 ! write (iout,*) 'chi(j,i)= ', chicat(j,1)
3307 ! write (iout,*) 'chip(i,j)= ', chippcat(1,j)
3308 ! write (iout,*) 'chip(j,i)= ', chippcat(j,1)
3309 ! write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
3310 ! write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
3311 ! write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
3312 ! write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
3313 ! write (iout,*) 'sig1= ', sigmap1cat(1,j)
3314 ! write (iout,*) 'chis(i,j)= ', chiscat(1,j)
3315 ! write (iout,*) 'chis(j,i)= ', chiscat(j,1)
3316 ! write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
3317 ! write (iout,*) 'a1= ', rborncat(j,1)
3318 ! write (iout,*) 'a2= ', rborncat(1,j)
3319 ! write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
3320 ! write (iout,*) 'alphapol1= ', alphapolcat(1,j)
3321 ! write (iout,*) 'w1= ', wqdipcat(1,1,j)
3322 ! write (iout,*) 'w2= ', wqdipcat(2,1,j)
3326 ! If ((i.eq.1).and.(j.eq.27)) then
3327 ! write (iout,*) 'SEP'
3328 ! Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3329 ! Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3334 allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
3336 do j=1,ntyp_molec(5)
3340 sigeps=dsign(1.0D0,epsij)
3342 aa_aq_cat(i,j)=epsij*rrij*rrij
3343 bb_aq_cat(i,j)=-sigeps*epsij*rrij
3348 do j=1,ntyp_molec(5)
3350 write (iout,*) 'i= ', i, ' j= ', j
3351 write (iout,*) 'epsi0= ', epscat(i,j)
3352 write (iout,*) 'sigma0= ', sigmacat(i,j)
3353 write (iout,*) 'chi1= ', chi1cat(i,j)
3354 write (iout,*) 'chi1= ', chi2cat(i,j)
3355 write (iout,*) 'chip1= ', chipp1cat(1,j)
3356 write (iout,*) 'chip2= ', chipp2cat(1,j)
3357 write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
3358 write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
3359 write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
3360 write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
3361 write (iout,*) 'sig1= ', sigmap1cat(1,j)
3362 write (iout,*) 'sig2= ', sigmap2cat(1,j)
3363 write (iout,*) 'chis1= ', chis1cat(1,j)
3364 write (iout,*) 'chis1= ', chis2cat(1,j)
3365 write (iout,*) 'nstatecat(i,j)= ', nstatecat(1,j)
3366 write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,1,j)
3367 write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
3368 write (iout,*) 'dhead2= ', dheadcat(1,2,1,j)
3369 write (iout,*) 'a1= ', rborn1cat(i,j)
3370 write (iout,*) 'a2= ', rborn2cat(i,j)
3371 write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
3372 write (iout,*) 'alphapol1= ', alphapolcat(1,j)
3373 write (iout,*) 'alphapol2= ', alphapolcat(j,1)
3374 write (iout,*) 'w1= ', wqdipcat(1,i,j)
3375 write (iout,*) 'w2= ', wqdipcat(2,i,j)
3376 write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(1,j)
3379 If ((i.eq.1).and.(j.eq.27)) then
3380 write (iout,*) 'SEP'
3381 Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3382 Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3392 write (iout,'(/a)') "Disulfide bridge parameters:"
3393 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3394 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3395 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3396 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3399 if (shield_mode.gt.0) then
3400 pi=4.0D0*datan(1.0D0)
3401 !C VSolvSphere the volume of solving sphere
3403 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3404 !C there will be no distinction between proline peptide group and normal peptide
3405 !C group in case of shielding parameters
3406 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3407 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3408 write (iout,*) VSolvSphere,VSolvSphere_div
3409 !C long axis of side chain
3411 long_r_sidechain(i)=vbldsc0(1,i)
3412 ! if (scelemode.eq.0) then
3413 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
3414 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
3416 ! short_r_sidechain(i)=sigma(i,i)
3418 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3425 111 write (iout,*) "Error reading bending energy parameters."
3427 112 write (iout,*) "Error reading rotamer energy parameters."
3429 113 write (iout,*) "Error reading torsional energy parameters."
3431 114 write (iout,*) "Error reading double torsional energy parameters."
3433 115 write (iout,*) &
3434 "Error reading cumulant (multibody energy) parameters."
3436 116 write (iout,*) "Error reading electrostatic energy parameters."
3438 117 write (iout,*) "Error reading side chain interaction parameters."
3440 118 write (iout,*) "Error reading SCp interaction parameters."
3442 119 write (iout,*) "Error reading SCCOR parameters"
3444 121 write (iout,*) "Error in Czybyshev parameters"
3447 call MPI_Finalize(Ierror)
3451 end subroutine parmread
3453 !-----------------------------------------------------------------------------
3455 !-----------------------------------------------------------------------------
3456 subroutine printmat(ldim,m,n,iout,key,a)
3459 character(len=3),dimension(n) :: key
3460 real(kind=8),dimension(ldim,n) :: a
3462 integer :: i,j,k,m,iout,nlim
3466 write (iout,1000) (key(k),k=i,nlim)
3468 1000 format (/5x,8(6x,a3))
3469 1020 format (/80(1h-)/)
3471 write (iout,1010) key(j),(a(j,k),k=i,nlim)
3474 1010 format (a3,2x,8(f9.4))
3476 end subroutine printmat
3477 !-----------------------------------------------------------------------------
3479 !-----------------------------------------------------------------------------
3481 ! Read the PDB file and convert the peptide geometry into virtual-chain
3484 use energy_data, only: itype,istype
3488 ! use control, only: rescode,sugarcode
3489 ! implicit real*8 (a-h,o-z)
3490 ! include 'DIMENSIONS'
3491 ! include 'COMMON.LOCAL'
3492 ! include 'COMMON.VAR'
3493 ! include 'COMMON.CHAIN'
3494 ! include 'COMMON.INTERACT'
3495 ! include 'COMMON.IOUNITS'
3496 ! include 'COMMON.GEO'
3497 ! include 'COMMON.NAMES'
3498 ! include 'COMMON.CONTROL'
3499 ! include 'COMMON.DISTFIT'
3500 ! include 'COMMON.SETUP'
3501 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
3503 logical :: lprn=.true.,fail
3504 real(kind=8),dimension(3) :: e1,e2,e3
3505 real(kind=8) :: dcj,efree_temp
3506 character(len=3) :: seq,res,res2
3507 character(len=5) :: atom
3508 character(len=80) :: card
3509 real(kind=8),dimension(3,20) :: sccor
3510 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
3511 integer :: isugar,molecprev,firstion
3512 character*1 :: sugar
3514 real(kind=8),dimension(3) :: ccc
3516 integer,dimension(2,maxres/3) :: hfrag_alloc
3517 integer,dimension(4,maxres/3) :: bfrag_alloc
3518 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
3519 real(kind=8),dimension(:,:), allocatable :: c_temporary
3520 integer,dimension(:,:) , allocatable :: itype_temporary
3521 integer,dimension(:),allocatable :: istype_temp
3528 ! write (2,*) "UNRES_PDB",unres_pdb
3548 !-----------------------------
3549 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
3550 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
3551 if(.not. allocated(istype)) allocate(istype(maxres))
3553 read (ipdbin,'(a80)',end=10) card
3554 write (iout,'(a)') card
3555 if (card(:5).eq.'HELIX') then
3558 read(card(22:25),*) hfrag(1,nhfrag)
3559 read(card(34:37),*) hfrag(2,nhfrag)
3561 if (card(:5).eq.'SHEET') then
3564 read(card(24:26),*) bfrag(1,nbfrag)
3565 read(card(35:37),*) bfrag(2,nbfrag)
3566 !rc----------------------------------------
3567 !rc to be corrected !!!
3568 bfrag(3,nbfrag)=bfrag(1,nbfrag)
3569 bfrag(4,nbfrag)=bfrag(2,nbfrag)
3570 !rc----------------------------------------
3572 if (card(:3).eq.'END') then
3574 else if (card(:3).eq.'TER') then
3579 itype(ires_old,molecule)=ntyp1_molec(molecule)
3580 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
3581 nres_molec(molecule)=nres_molec(molecule)+2
3583 ! write (iout,*) "Chain ended",ires,ishift,ires_old
3586 dc(j,ires)=sccor(j,iii)
3589 call sccenter(ires,iii,sccor)
3595 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
3596 ! Fish out the ATOM cards.
3597 ! write(iout,*) 'card',card(1:20)
3598 ! print *,"ATU ",card(1:6), CARD(3:6)
3600 if (index(card(1:4),'ATOM').gt.0) then
3601 read (card(12:16),*) atom
3602 ! write (iout,*) "! ",atom," !",ires
3603 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
3604 read (card(23:26),*) ires
3605 read (card(18:20),'(a3)') res
3606 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
3607 ! & " ires_old",ires_old
3608 ! write (iout,*) "ishift",ishift," ishift1",ishift1
3609 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
3610 if (ires-ishift+ishift1.ne.ires_old) then
3611 ! Calculate the CM of the preceding residue.
3612 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
3614 ! write (iout,*) "Calculating sidechain center iii",iii
3617 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
3620 call sccenter(ires_old,iii,sccor)
3624 ! Start new residue.
3625 if (res.eq.'Cl-' .or. res.eq.'Na+') then
3628 else if (ibeg.eq.1) then
3629 write (iout,*) "BEG ires",ires
3631 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3634 nres_molec(molecule)=nres_molec(molecule)+1
3636 ires=ires-ishift+ishift1
3638 ! write (iout,*) "ishift",ishift," ires",ires,&
3639 ! " ires_old",ires_old
3641 else if (ibeg.eq.2) then
3643 ishift=-ires_old+ires-1 !!!!!
3644 ishift1=ishift1-1 !!!!!
3645 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3646 ires=ires-ishift+ishift1
3647 ! print *,ires,ishift,ishift1
3651 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
3652 ires=ires-ishift+ishift1
3655 ! print *,'atom',ires,atom
3656 if (res.eq.'ACE' .or. res.eq.'NHE') then
3659 if (atom.eq.'CA '.or.atom.eq.'N ') then
3661 itype(ires,molecule)=rescode(ires,res,0,molecule)
3663 ! nres_molec(molecule)=nres_molec(molecule)+1
3667 itype(ires,molecule)=rescode(ires,res2,0,molecule)
3668 ! nres_molec(molecule)=nres_molec(molecule)+1
3669 read (card(19:19),'(a1)') sugar
3670 isugar=sugarcode(sugar,ires)
3671 ! if (ibeg.eq.1) then
3675 ! print *,"ires=",ires,istype(ires)
3681 ires=ires-ishift+ishift1
3683 ! write (iout,*) "ires_old",ires_old," ires",ires
3684 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3687 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3688 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
3689 res.eq.'NHE'.and.atom(:2).eq.'HN') then
3690 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3691 ! print *,ires,ishift,ishift1
3692 ! write (iout,*) "backbone ",atom
3694 write (iout,'(2i3,2x,a,3f8.3)') &
3695 ires,itype(ires,1),res,(c(j,ires),j=1,3)
3698 nres_molec(molecule)=nres_molec(molecule)+1
3700 sccor(j,iii)=c(j,ires)
3702 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
3703 atom.eq."C2'" .or. atom.eq."C3'" &
3704 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
3705 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
3706 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
3707 ! print *,ires,ishift,ishift1
3711 ! sccor(j,iii)=c(j,ires)
3714 c(j,ires)=c(j,ires)+ccc(j)/5.0
3716 print *,counter,molecule
3717 if (counter.eq.5) then
3719 nres_molec(molecule)=nres_molec(molecule)+1
3722 ! sccor(j,iii)=c(j,ires)
3726 ! print *, "ATOM",atom(1:3)
3727 else if (atom.eq."C5'") then
3728 read (card(19:19),'(a1)') sugar
3729 isugar=sugarcode(sugar,ires)
3734 ! print *,ires,istype(ires)
3738 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
3739 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3740 nres_molec(molecule)=nres_molec(molecule)+1
3741 print *,"nres_molec(molecule)",nres_molec(molecule),ires
3745 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3747 else if ((atom.eq."C1'").and.unres_pdb) then
3749 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3750 ! write (*,*) card(23:27),ires,itype(ires,1)
3751 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
3752 atom.ne.'N' .and. atom.ne.'C' .and. &
3753 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
3754 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
3755 .and. atom.ne.'P '.and. &
3756 atom(1:1).ne.'H' .and. &
3757 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
3759 ! write (iout,*) "sidechain ",atom
3760 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3761 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
3762 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3764 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3767 ! print *,"IONS",ions,card(1:6)
3768 else if ((ions).and.(card(1:6).eq.'HETATM')) then
3769 if (firstion.eq.0) then
3773 dc(j,ires)=sccor(j,iii)
3776 call sccenter(ires,iii,sccor)
3779 read (card(12:16),*) atom
3780 ! print *,"HETATOM", atom
3781 read (card(18:20),'(a3)') res
3782 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3783 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3784 .or.(atom(1:2).eq.'K ')) &
3787 if (molecule.ne.5) molecprev=molecule
3789 nres_molec(molecule)=nres_molec(molecule)+1
3790 print *,"HERE",nres_molec(molecule)
3792 itype(ires,molecule)=rescode(ires,res,0,molecule)
3793 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3797 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3798 if (ires.eq.0) return
3799 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3802 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
3804 nres_molec(molecule)=nres_molec(molecule)-2
3805 print *,'I have',nres, nres_molec(:)
3807 do k=1,4 ! ions are without dummy
3808 if (nres_molec(k).eq.0) cycle
3810 ! write (iout,*) i,itype(i,1)
3811 ! if (itype(i,1).eq.ntyp1) then
3812 ! write (iout,*) "dummy",i,itype(i,1)
3814 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3815 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3819 if (itype(i,k).eq.ntyp1_molec(k)) then
3820 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3821 if (itype(i+2,k).eq.0) then
3822 ! print *,"masz sieczke"
3824 if (itype(i+2,j).ne.0) then
3826 itype(i+1,j)=ntyp1_molec(j)
3827 nres_molec(k)=nres_molec(k)-1
3828 nres_molec(j)=nres_molec(j)+1
3834 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3835 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3836 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3837 ! if (unres_pdb) then
3838 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3839 ! print *,i,'tu dochodze'
3840 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3848 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3852 dcj=(c(j,i-2)-c(j,i-3))/2.0
3853 if (dcj.eq.0) dcj=1.23591524223
3858 else !itype(i+1,1).eq.ntyp1
3859 ! if (unres_pdb) then
3860 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3861 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3868 ! c(j,i)=c(j,i+1)-1.9d0*e2(j)
3869 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
3873 dcj=(c(j,i+3)-c(j,i+2))/2.0
3874 if (dcj.eq.0) dcj=1.23591524223
3879 endif !itype(i+1,1).eq.ntyp1
3880 endif !itype.eq.ntyp1
3884 ! Calculate the CM of the last side chain.
3888 dc(j,ires)=sccor(j,iii)
3891 call sccenter(ires,iii,sccor)
3897 ! print *,"molecule",molecule
3898 if ((itype(nres,1).ne.10)) then
3900 if (molecule.eq.5) molecule=molecprev
3901 itype(nres,molecule)=ntyp1_molec(molecule)
3902 nres_molec(molecule)=nres_molec(molecule)+1
3903 ! if (unres_pdb) then
3904 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3905 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3912 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3916 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3917 c(j,nres)=c(j,nres-1)+dcj
3918 c(j,2*nres)=c(j,nres)
3922 ! print *,'I have',nres, nres_molec(:)
3924 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3926 if (nres.ne.nres0) then
3927 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3929 stop "Error nres value in WHAM input"
3932 !---------------------------------
3933 !el reallocate tables
3936 ! hfrag_alloc(j,i)=hfrag(j,i)
3939 ! bfrag_alloc(j,i)=bfrag(j,i)
3945 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3946 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3947 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3948 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3952 ! hfrag(j,i)=hfrag_alloc(j,i)
3957 ! bfrag(j,i)=bfrag_alloc(j,i)
3960 !el end reallocate tables
3961 !---------------------------------
3969 c(j,2*nres)=c(j,nres)
3972 if (itype(1,1).eq.ntyp1) then
3976 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3977 call refsys(2,3,4,e1,e2,e3,fail)
3984 ! c(j,1)=c(j,2)-1.9d0*e2(j)
3985 c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0)
3989 dcj=(c(j,4)-c(j,3))/2.0
3995 ! First lets assign correct dummy to correct type of chain
3997 if (itype(1,1).eq.ntyp1) then
3998 if (itype(2,1).eq.0) then
4000 if (itype(2,j).ne.0) then
4002 itype(1,j)=ntyp1_molec(j)
4003 nres_molec(1)=nres_molec(1)-1
4004 nres_molec(j)=nres_molec(j)+1
4011 print *,'I have',nres, nres_molec(:)
4013 ! Copy the coordinates to reference coordinates
4019 ! Calculate internal coordinates.
4021 write (iout,'(/a)') &
4022 "Cartesian coordinates of the reference structure"
4023 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4024 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4026 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4027 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4028 (c(j,ires+nres),j=1,3)
4031 ! znamy już nres wiec mozna alokowac tablice
4032 ! Calculate internal coordinates.
4033 if(me.eq.king.or..not.out1file)then
4034 write (iout,'(a)') &
4035 "Backbone and SC coordinates as read from the PDB"
4037 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
4038 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
4039 (c(j,nres+ires),j=1,3)
4042 ! NOW LETS ROCK! SORTING
4043 allocate(c_temporary(3,2*nres))
4044 allocate(itype_temporary(nres,5))
4045 if (.not.allocated(molnum)) allocate(molnum(nres+1))
4046 if (.not.allocated(istype)) write(iout,*) &
4047 "SOMETHING WRONG WITH ISTYTPE"
4048 allocate(istype_temp(nres))
4049 itype_temporary(:,:)=0
4053 if (itype(i,k).ne.0) then
4055 c_temporary(j,seqalingbegin)=c(j,i)
4056 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
4059 itype_temporary(seqalingbegin,k)=itype(i,k)
4060 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
4061 istype_temp(seqalingbegin)=istype(i)
4062 molnum(seqalingbegin)=k
4063 seqalingbegin=seqalingbegin+1
4069 c(j,i)=c_temporary(j,i)
4074 itype(i,k)=itype_temporary(i,k)
4075 istype(i)=istype_temp(i)
4078 if ((itype(1,1).eq.ntyp1).and.itype(2,5).ne.0) then
4079 ! I have only ions now dummy atoms in the system
4081 itype(1,5)=itype(2,5)
4084 itype(i,5)=itype(i+1,5)
4088 nres_molec(1)=nres_molec(1)-1
4090 ! if (itype(1,1).eq.ntyp1) then
4093 ! if (unres_pdb) then
4094 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4095 ! call refsys(2,3,4,e1,e2,e3,fail)
4102 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4106 ! dcj=(c(j,4)-c(j,3))/2.0
4108 ! c(j,nres+1)=c(j,1)
4114 write (iout,'(/a)') &
4115 "Cartesian coordinates of the reference structure after sorting"
4116 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4117 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4119 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4120 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4121 (c(j,ires+nres),j=1,3)
4125 ! print *,seqalingbegin,nres
4126 if(.not.allocated(vbld)) then
4127 allocate(vbld(2*nres))
4132 if(.not.allocated(vbld_inv)) then
4133 allocate(vbld_inv(2*nres))
4139 if(.not.allocated(theta)) then
4140 allocate(theta(nres+2))
4144 if(.not.allocated(phi)) allocate(phi(nres+2))
4145 if(.not.allocated(alph)) allocate(alph(nres+2))
4146 if(.not.allocated(omeg)) allocate(omeg(nres+2))
4147 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
4148 if(.not.allocated(phiref)) allocate(phiref(nres+2))
4149 if(.not.allocated(costtab)) allocate(costtab(nres))
4150 if(.not.allocated(sinttab)) allocate(sinttab(nres))
4151 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
4152 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
4153 if(.not.allocated(xxref)) allocate(xxref(nres))
4154 if(.not.allocated(yyref)) allocate(yyref(nres))
4155 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
4156 if(.not.allocated(dc_norm)) then
4157 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
4158 allocate(dc_norm(3,0:2*nres+2))
4162 call int_from_cart(.true.,.false.)
4163 call sc_loc_geom(.false.)
4165 thetaref(i)=theta(i)
4175 dc(j,i)=c(j,i+1)-c(j,i)
4176 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
4181 dc(j,i+nres)=c(j,i+nres)-c(j,i)
4182 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
4184 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
4188 ! Copy the coordinates to reference coordinates
4189 ! Splits to single chain if occurs
4193 ! cref(j,i,cou)=c(j,i)
4197 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
4198 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
4199 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
4200 !-----------------------------
4204 write (iout,*) "symetr", symetr
4207 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4209 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
4212 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4217 cref(j,i,cou)=c(j,i)
4218 cref(j,i+nres,cou)=c(j,i+nres)
4220 chain_rep(j,lll,kkk)=c(j,i)
4221 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
4225 write (iout,*) chain_length
4226 if (chain_length.eq.0) chain_length=nres
4228 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
4229 chain_rep(j,chain_length+nres,symetr) &
4230 =chain_rep(j,chain_length+nres,1)
4233 ! write (iout,*) "spraw lancuchy",chain_length,symetr
4235 ! do kkk=1,chain_length
4236 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
4240 ! makes copy of chains
4241 write (iout,*) "symetr", symetr
4246 if (symetr.gt.1) then
4253 write(iout,*) (tabperm(i,kkk),kkk=1,4)
4259 write (iout,*) i,icha
4260 do lll=1,chain_length
4262 if (cou.le.nres) then
4264 kupa=mod(lll,chain_length)
4265 iprzes=(kkk-1)*chain_length+lll
4266 if (kupa.eq.0) kupa=chain_length
4267 write (iout,*) "kupa", kupa
4268 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
4269 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
4276 !-koniec robienia kopii
4279 write (iout,*) "nowa struktura", nperm
4281 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
4283 cref(3,i,kkk),cref(1,nres+i,kkk),&
4284 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
4286 100 format (//' alpha-carbon coordinates ',&
4287 ' centroid coordinates'/ &
4288 ' ', 6X,'X',11X,'Y',11X,'Z', &
4289 10X,'X',11X,'Y',11X,'Z')
4290 110 format (a,'(',i5,')',6f12.5)
4296 bfrag(i,j)=bfrag(i,j)-ishift
4302 hfrag(i,j)=hfrag(i,j)-ishift
4308 end subroutine readpdb
4309 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
4310 !-----------------------------------------------------------------------------
4312 !-----------------------------------------------------------------------------
4313 subroutine read_control
4327 use random, only: random_init
4328 ! implicit real*8 (a-h,o-z)
4329 ! include 'DIMENSIONS'
4331 use prng, only:prng_restart
4333 logical :: OKRandom!, prng_restart
4336 ! include 'COMMON.IOUNITS'
4337 ! include 'COMMON.TIME1'
4338 ! include 'COMMON.THREAD'
4339 ! include 'COMMON.SBRIDGE'
4340 ! include 'COMMON.CONTROL'
4341 ! include 'COMMON.MCM'
4342 ! include 'COMMON.MAP'
4343 ! include 'COMMON.HEADER'
4344 ! include 'COMMON.CSA'
4345 ! include 'COMMON.CHAIN'
4346 ! include 'COMMON.MUCA'
4347 ! include 'COMMON.MD'
4348 ! include 'COMMON.FFIELD'
4349 ! include 'COMMON.INTERACT'
4350 ! include 'COMMON.SETUP'
4351 !el integer :: KDIAG,ICORFL,IXDR
4352 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
4353 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
4354 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
4355 ! character(len=80) :: ucase
4356 character(len=640) :: controlcard
4358 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
4364 read (INP,'(a)') titel
4365 call card_concat(controlcard,.true.)
4366 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
4367 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
4368 call reada(controlcard,'SEED',seed,0.0D0)
4369 call random_init(seed)
4370 ! Set up the time limit (caution! The time must be input in minutes!)
4371 read_cart=index(controlcard,'READ_CART').gt.0
4372 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
4373 call readi(controlcard,'SYM',symetr,1)
4374 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
4375 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
4376 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
4377 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
4378 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
4379 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
4380 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
4381 call reada(controlcard,'DRMS',drms,0.1D0)
4382 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4383 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
4384 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
4385 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
4386 write (iout,'(a,f10.1)')'DRMS = ',drms
4387 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
4388 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
4390 call readi(controlcard,'NZ_START',nz_start,0)
4391 call readi(controlcard,'NZ_END',nz_end,0)
4392 call readi(controlcard,'IZ_SC',iz_sc,0)
4393 timlim=60.0D0*timlim
4394 safety = 60.0d0*safety
4397 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4398 !C SHIELD keyword sets if the shielding effect of side-chains is used
4399 !C 0 denots no shielding is used all peptide are equally despite the
4400 !C solvent accesible area
4401 !C 1 the newly introduced function
4402 !C 2 reseved for further possible developement
4403 call readi(controlcard,'SHIELD',shield_mode,0)
4404 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4405 write(iout,*) "shield_mode",shield_mode
4406 call readi(controlcard,'TORMODE',tor_mode,0)
4407 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4408 write(iout,*) "torsional and valence angle mode",tor_mode
4410 !C Varibles set size of box
4411 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
4412 protein=index(controlcard,"PROTEIN").gt.0
4413 ions=index(controlcard,"IONS").gt.0
4414 call readi(controlcard,'OLDION',oldion,1)
4415 nucleic=index(controlcard,"NUCLEIC").gt.0
4416 write (iout,*) "with_theta_constr ",with_theta_constr
4417 AFMlog=(index(controlcard,'AFM'))
4418 selfguide=(index(controlcard,'SELFGUIDE'))
4419 print *,'AFMlog',AFMlog,selfguide,"KUPA"
4420 call readi(controlcard,'GENCONSTR',genconstr,0)
4421 call reada(controlcard,'BOXX',boxxsize,100.0d0)
4422 call reada(controlcard,'BOXY',boxysize,100.0d0)
4423 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
4424 call readi(controlcard,'TUBEMOD',tubemode,0)
4425 print *,"SCELE",scelemode
4426 call readi(controlcard,"SCELEMODE",scelemode,0)
4427 print *,"SCELE",scelemode
4429 ! elemode = 0 is orignal UNRES electrostatics
4430 ! elemode = 1 is "Momo" potentials in progress
4431 ! elemode = 2 is in development EVALD
4434 write (iout,*) TUBEmode,"TUBEMODE"
4435 if (TUBEmode.gt.0) then
4436 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
4437 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
4438 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
4439 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
4440 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
4441 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
4442 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
4443 buftubebot=bordtubebot+tubebufthick
4444 buftubetop=bordtubetop-tubebufthick
4447 ! CUTOFFF ON ELECTROSTATICS
4448 call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
4449 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
4450 write(iout,*) "R_CUT_ELE=",r_cut_ele
4451 ! Lipidic parameters
4452 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
4453 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
4454 if (lipthick.gt.0.0d0) then
4455 bordliptop=(boxzsize+lipthick)/2.0
4456 bordlipbot=bordliptop-lipthick
4457 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
4458 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
4459 buflipbot=bordlipbot+lipbufthick
4460 bufliptop=bordliptop-lipbufthick
4461 if ((lipbufthick*2.0d0).gt.lipthick) &
4462 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
4463 endif !lipthick.gt.0
4464 write(iout,*) "bordliptop=",bordliptop
4465 write(iout,*) "bordlipbot=",bordlipbot
4466 write(iout,*) "bufliptop=",bufliptop
4467 write(iout,*) "buflipbot=",buflipbot
4468 write (iout,*) "SHIELD MODE",shield_mode
4470 !C-------------------------
4471 minim=(index(controlcard,'MINIMIZE').gt.0)
4472 dccart=(index(controlcard,'CART').gt.0)
4473 overlapsc=(index(controlcard,'OVERLAP').gt.0)
4474 overlapsc=.not.overlapsc
4475 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
4476 searchsc=.not.searchsc
4477 sideadd=(index(controlcard,'SIDEADD').gt.0)
4478 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
4479 outpdb=(index(controlcard,'PDBOUT').gt.0)
4480 outmol2=(index(controlcard,'MOL2OUT').gt.0)
4481 pdbref=(index(controlcard,'PDBREF').gt.0)
4482 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
4483 indpdb=index(controlcard,'PDBSTART')
4484 extconf=(index(controlcard,'EXTCONF').gt.0)
4485 call readi(controlcard,'IPRINT',iprint,0)
4486 call readi(controlcard,'MAXGEN',maxgen,10000)
4487 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
4488 call readi(controlcard,"KDIAG",kdiag,0)
4489 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
4490 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
4491 write (iout,*) "RESCALE_MODE",rescale_mode
4492 split_ene=index(controlcard,'SPLIT_ENE').gt.0
4493 if (index(controlcard,'REGULAR').gt.0.0D0) then
4494 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4498 if (index(controlcard,'CHECKGRAD').gt.0) then
4500 if (index(controlcard,'CART').gt.0) then
4502 elseif (index(controlcard,'CARINT').gt.0) then
4507 elseif (index(controlcard,'THREAD').gt.0) then
4509 call readi(controlcard,'THREAD',nthread,0)
4510 if (nthread.gt.0) then
4511 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4514 write (iout,'(a)')'A number has to follow the THREAD keyword.'
4515 stop 'Error termination in Read_Control.'
4517 else if (index(controlcard,'MCMA').gt.0) then
4519 else if (index(controlcard,'MCEE').gt.0) then
4521 else if (index(controlcard,'MULTCONF').gt.0) then
4523 else if (index(controlcard,'MAP').gt.0) then
4525 call readi(controlcard,'MAP',nmap,0)
4526 else if (index(controlcard,'CSA').gt.0) then
4528 !rc else if (index(controlcard,'ZSCORE').gt.0) then
4530 !rc ZSCORE is rm from UNRES, modecalc=9 is available
4533 !fcm else if (index(controlcard,'MCMF').gt.0) then
4535 else if (index(controlcard,'SOFTREG').gt.0) then
4537 else if (index(controlcard,'CHECK_BOND').gt.0) then
4539 else if (index(controlcard,'TEST').gt.0) then
4541 else if (index(controlcard,'MD').gt.0) then
4543 else if (index(controlcard,'RE ').gt.0) then
4547 lmuca=index(controlcard,'MUCA').gt.0
4548 call readi(controlcard,'MUCADYN',mucadyn,0)
4549 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
4550 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
4552 write (iout,*) 'MUCADYN=',mucadyn
4553 write (iout,*) 'MUCASMOOTH=',muca_smooth
4556 iscode=index(controlcard,'ONE_LETTER')
4557 indphi=index(controlcard,'PHI')
4558 indback=index(controlcard,'BACK')
4559 iranconf=index(controlcard,'RAND_CONF')
4560 i2ndstr=index(controlcard,'USE_SEC_PRED')
4561 gradout=index(controlcard,'GRADOUT').gt.0
4562 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
4563 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
4564 if (me.eq.king .or. .not.out1file ) &
4565 write (iout,*) "DISTCHAINMAX",distchainmax
4567 if(me.eq.king.or..not.out1file) &
4568 write (iout,'(2a)') diagmeth(kdiag),&
4569 ' routine used to diagonalize matrices.'
4570 if (shield_mode.gt.0) then
4571 pi=4.0D0*datan(1.0D0)
4572 !C VSolvSphere the volume of solving sphere
4574 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
4575 !C there will be no distinction between proline peptide group and normal peptide
4576 !C group in case of shielding parameters
4577 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
4578 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
4579 write (iout,*) VSolvSphere,VSolvSphere_div
4580 !C long axis of side chain
4582 ! long_r_sidechain(i)=vbldsc0(1,i)
4583 ! short_r_sidechain(i)=sigma0(i)
4584 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
4591 end subroutine read_control
4592 !-----------------------------------------------------------------------------
4593 subroutine read_REMDpar
4595 ! Read REMD settings
4602 use control_data, only:out1file
4603 ! implicit real*8 (a-h,o-z)
4604 ! include 'DIMENSIONS'
4605 ! include 'COMMON.IOUNITS'
4606 ! include 'COMMON.TIME1'
4607 ! include 'COMMON.MD'
4610 !el include 'COMMON.LANGEVIN'
4612 !el include 'COMMON.LANGEVIN.lang0'
4614 ! include 'COMMON.INTERACT'
4615 ! include 'COMMON.NAMES'
4616 ! include 'COMMON.GEO'
4617 ! include 'COMMON.REMD'
4618 ! include 'COMMON.CONTROL'
4619 ! include 'COMMON.SETUP'
4620 ! character(len=80) :: ucase
4621 character(len=320) :: controlcard
4622 character(len=3200) :: controlcard1
4623 integer :: iremd_m_total
4626 ! real(kind=8) :: var,ene
4628 if(me.eq.king.or..not.out1file) &
4629 write (iout,*) "REMD setup"
4631 call card_concat(controlcard,.true.)
4632 call readi(controlcard,"NREP",nrep,3)
4633 call readi(controlcard,"NSTEX",nstex,1000)
4634 call reada(controlcard,"RETMIN",retmin,10.0d0)
4635 call reada(controlcard,"RETMAX",retmax,1000.0d0)
4636 mremdsync=(index(controlcard,'SYNC').gt.0)
4637 call readi(controlcard,"NSYN",i_sync_step,100)
4638 restart1file=(index(controlcard,'REST1FILE').gt.0)
4639 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
4640 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
4641 if(max_cache_traj_use.gt.max_cache_traj) &
4642 max_cache_traj_use=max_cache_traj
4643 if(me.eq.king.or..not.out1file) then
4644 !d if (traj1file) then
4645 !rc caching is in testing - NTWX is not ignored
4646 !d write (iout,*) "NTWX value is ignored"
4647 !d write (iout,*) " trajectory is stored to one file by master"
4648 !d write (iout,*) " before exchange at NSTEX intervals"
4650 write (iout,*) "NREP= ",nrep
4651 write (iout,*) "NSTEX= ",nstex
4652 write (iout,*) "SYNC= ",mremdsync
4653 write (iout,*) "NSYN= ",i_sync_step
4654 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
4657 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
4658 if (index(controlcard,'TLIST').gt.0) then
4660 call card_concat(controlcard1,.true.)
4661 read(controlcard1,*) (remd_t(i),i=1,nrep)
4662 if(me.eq.king.or..not.out1file) &
4663 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
4666 if (index(controlcard,'MLIST').gt.0) then
4668 call card_concat(controlcard1,.true.)
4669 read(controlcard1,*) (remd_m(i),i=1,nrep)
4670 if(me.eq.king.or..not.out1file) then
4671 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
4674 iremd_m_total=iremd_m_total+remd_m(i)
4676 write (iout,*) 'Total number of replicas ',iremd_m_total
4679 if(me.eq.king.or..not.out1file) &
4680 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
4682 end subroutine read_REMDpar
4683 !-----------------------------------------------------------------------------
4684 subroutine read_MDpar
4688 use control_data, only: r_cut,rlamb,out1file
4690 use geometry_data, only: pi
4692 ! implicit real*8 (a-h,o-z)
4693 ! include 'DIMENSIONS'
4694 ! include 'COMMON.IOUNITS'
4695 ! include 'COMMON.TIME1'
4696 ! include 'COMMON.MD'
4699 !el include 'COMMON.LANGEVIN'
4701 !el include 'COMMON.LANGEVIN.lang0'
4703 ! include 'COMMON.INTERACT'
4704 ! include 'COMMON.NAMES'
4705 ! include 'COMMON.GEO'
4706 ! include 'COMMON.SETUP'
4707 ! include 'COMMON.CONTROL'
4708 ! include 'COMMON.SPLITELE'
4709 ! character(len=80) :: ucase
4710 character(len=320) :: controlcard
4715 call card_concat(controlcard,.true.)
4716 call readi(controlcard,"NSTEP",n_timestep,1000000)
4717 call readi(controlcard,"NTWE",ntwe,100)
4718 call readi(controlcard,"NTWX",ntwx,1000)
4719 call reada(controlcard,"DT",d_time,1.0d-1)
4720 call reada(controlcard,"DVMAX",dvmax,2.0d1)
4721 call reada(controlcard,"DAMAX",damax,1.0d1)
4722 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
4723 call readi(controlcard,"LANG",lang,0)
4724 RESPA = index(controlcard,"RESPA") .gt. 0
4725 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
4726 ntime_split0=ntime_split
4727 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
4728 ntime_split0=ntime_split
4729 call reada(controlcard,"R_CUT",r_cut,2.0d0)
4730 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
4731 rest = index(controlcard,"REST").gt.0
4732 tbf = index(controlcard,"TBF").gt.0
4733 usampl = index(controlcard,"USAMPL").gt.0
4734 mdpdb = index(controlcard,"MDPDB").gt.0
4735 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4736 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
4737 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
4738 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
4739 if (count_reset_moment.eq.0) count_reset_moment=1000000000
4740 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
4741 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
4742 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
4743 if (count_reset_vel.eq.0) count_reset_vel=1000000000
4744 large = index(controlcard,"LARGE").gt.0
4745 print_compon = index(controlcard,"PRINT_COMPON").gt.0
4746 rattle = index(controlcard,"RATTLE").gt.0
4747 preminim=(index(controlcard,'PREMINIM').gt.0)
4748 write (iout,*) "PREMINIM ",preminim
4749 dccart=(index(controlcard,'CART').gt.0)
4750 if (preminim) call read_minim
4751 ! if performing umbrella sampling, fragments constrained are read from the fragment file
4757 if(me.eq.king.or..not.out1file) then
4759 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
4761 write (iout,'(a)') "The units are:"
4762 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
4763 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
4764 " acceleration: angstrom/(48.9 fs)**2"
4765 write (iout,'(a)') "energy: kcal/mol, temperature: K"
4767 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
4768 write (iout,'(a60,f10.5,a)') &
4769 "Initial time step of numerical integration:",d_time,&
4771 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
4773 write (iout,'(2a,i4,a)') &
4774 "A-MTS algorithm used; initial time step for fast-varying",&
4775 " short-range forces split into",ntime_split," steps."
4776 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
4777 r_cut," lambda",rlamb
4779 write (iout,'(2a,f10.5)') &
4780 "Maximum acceleration threshold to reduce the time step",&
4781 "/increase split number:",damax
4782 write (iout,'(2a,f10.5)') &
4783 "Maximum predicted energy drift to reduce the timestep",&
4784 "/increase split number:",edriftmax
4785 write (iout,'(a60,f10.5)') &
4786 "Maximum velocity threshold to reduce velocities:",dvmax
4787 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
4788 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
4789 if (rattle) write (iout,'(a60)') &
4790 "Rattle algorithm used to constrain the virtual bonds"
4794 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
4795 call reada(controlcard,"RWAT",rwat,1.4d0)
4796 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
4797 surfarea=index(controlcard,"SURFAREA").gt.0
4798 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
4799 if(me.eq.king.or..not.out1file)then
4800 write (iout,'(/a,$)') "Langevin dynamics calculation"
4802 write (iout,'(a/)') &
4803 " with direct integration of Langevin equations"
4804 else if (lang.eq.2) then
4805 write (iout,'(a/)') " with TINKER stochasic MD integrator"
4806 else if (lang.eq.3) then
4807 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
4808 else if (lang.eq.4) then
4809 write (iout,'(a/)') " in overdamped mode"
4811 write (iout,'(//a,i5)') &
4812 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4815 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4816 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4817 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4818 write (iout,'(a60,f10.5)') &
4819 "Scaling factor of the friction forces:",scal_fric
4820 if (surfarea) write (iout,'(2a,i10,a)') &
4821 "Friction coefficients will be scaled by solvent-accessible",&
4822 " surface area every",reset_fricmat," steps."
4824 ! Calculate friction coefficients and bounds of stochastic forces
4825 eta=6*pi*cPoise*etawat
4826 if(me.eq.king.or..not.out1file) &
4827 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4830 do j=1,5 !types of molecules
4831 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4832 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4834 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4835 do j=1,5 !types of molecules
4837 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4838 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4842 if(me.eq.king.or..not.out1file)then
4843 write (iout,'(/2a/)') &
4844 "Radii of site types and friction coefficients and std's of",&
4845 " stochastic forces of fully exposed sites"
4846 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4848 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4849 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4853 if(me.eq.king.or..not.out1file)then
4854 write (iout,'(a)') "Berendsen bath calculation"
4855 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4856 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4858 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4859 count_reset_moment," steps"
4861 write (iout,'(a,i10,a)') &
4862 "Velocities will be reset at random every",count_reset_vel,&
4866 if(me.eq.king.or..not.out1file) &
4867 write (iout,'(a31)') "Microcanonical mode calculation"
4869 if(me.eq.king.or..not.out1file)then
4870 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4872 write(iout,*) "MD running with constraints."
4873 write(iout,*) "Equilibration time ", eq_time, " mtus."
4874 write(iout,*) "Constraining ", nfrag," fragments."
4875 write(iout,*) "Length of each fragment, weight and q0:"
4877 write (iout,*) "Set of restraints #",iset
4879 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4880 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4882 write(iout,*) "constraints between ", npair, "fragments."
4883 write(iout,*) "constraint pairs, weights and q0:"
4885 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4886 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4888 write(iout,*) "angle constraints within ", nfrag_back,&
4889 "backbone fragments."
4890 write(iout,*) "fragment, weights:"
4892 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4893 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4894 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4897 iset=mod(kolor,nset)+1
4900 if(me.eq.king.or..not.out1file) &
4901 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4903 end subroutine read_MDpar
4904 !-----------------------------------------------------------------------------
4908 ! implicit real*8 (a-h,o-z)
4909 ! include 'DIMENSIONS'
4910 ! include 'COMMON.MAP'
4911 ! include 'COMMON.IOUNITS'
4912 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4913 character(len=80) :: mapcard !,ucase
4916 ! real(kind=8) :: var,ene
4919 read (inp,'(a)') mapcard
4920 mapcard=ucase(mapcard)
4921 if (index(mapcard,'PHI').gt.0) then
4923 else if (index(mapcard,'THE').gt.0) then
4925 else if (index(mapcard,'ALP').gt.0) then
4927 else if (index(mapcard,'OME').gt.0) then
4930 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4931 stop 'Error - illegal variable spec in MAP card.'
4933 call readi (mapcard,'RES1',res1(imap),0)
4934 call readi (mapcard,'RES2',res2(imap),0)
4935 if (res1(imap).eq.0) then
4936 res1(imap)=res2(imap)
4937 else if (res2(imap).eq.0) then
4938 res2(imap)=res1(imap)
4940 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4941 write (iout,'(a)') &
4942 'Error - illegal definition of variable group in MAP.'
4943 stop 'Error - illegal definition of variable group in MAP.'
4945 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4946 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4947 call readi(mapcard,'NSTEP',nstep(imap),0)
4948 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4949 write (iout,'(a)') &
4950 'Illegal boundary and/or step size specification in MAP.'
4951 stop 'Illegal boundary and/or step size specification in MAP.'
4955 end subroutine map_read
4956 !-----------------------------------------------------------------------------
4959 use control_data, only: vdisulf
4961 ! implicit real*8 (a-h,o-z)
4962 ! include 'DIMENSIONS'
4963 ! include 'COMMON.IOUNITS'
4964 ! include 'COMMON.GEO'
4965 ! include 'COMMON.CSA'
4966 ! include 'COMMON.BANK'
4967 ! include 'COMMON.CONTROL'
4968 ! character(len=80) :: ucase
4969 character(len=620) :: mcmcard
4971 ! integer :: ntf,ik,iw_pdb
4972 ! real(kind=8) :: var,ene
4974 call card_concat(mcmcard,.true.)
4976 call readi(mcmcard,'NCONF',nconf,50)
4977 call readi(mcmcard,'NADD',nadd,0)
4978 call readi(mcmcard,'JSTART',jstart,1)
4979 call readi(mcmcard,'JEND',jend,1)
4980 call readi(mcmcard,'NSTMAX',nstmax,500000)
4981 call readi(mcmcard,'N0',n0,1)
4982 call readi(mcmcard,'N1',n1,6)
4983 call readi(mcmcard,'N2',n2,4)
4984 call readi(mcmcard,'N3',n3,0)
4985 call readi(mcmcard,'N4',n4,0)
4986 call readi(mcmcard,'N5',n5,0)
4987 call readi(mcmcard,'N6',n6,10)
4988 call readi(mcmcard,'N7',n7,0)
4989 call readi(mcmcard,'N8',n8,0)
4990 call readi(mcmcard,'N9',n9,0)
4991 call readi(mcmcard,'N14',n14,0)
4992 call readi(mcmcard,'N15',n15,0)
4993 call readi(mcmcard,'N16',n16,0)
4994 call readi(mcmcard,'N17',n17,0)
4995 call readi(mcmcard,'N18',n18,0)
4997 vdisulf=(index(mcmcard,'DYNSS').gt.0)
4999 call readi(mcmcard,'NDIFF',ndiff,2)
5000 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
5001 call readi(mcmcard,'IS1',is1,1)
5002 call readi(mcmcard,'IS2',is2,8)
5003 call readi(mcmcard,'NRAN0',nran0,4)
5004 call readi(mcmcard,'NRAN1',nran1,2)
5005 call readi(mcmcard,'IRR',irr,1)
5006 call readi(mcmcard,'NSEED',nseed,20)
5007 call readi(mcmcard,'NTOTAL',ntotal,10000)
5008 call reada(mcmcard,'CUT1',cut1,2.0d0)
5009 call reada(mcmcard,'CUT2',cut2,5.0d0)
5010 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
5011 call readi(mcmcard,'ICMAX',icmax,3)
5012 call readi(mcmcard,'IRESTART',irestart,0)
5013 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
5016 call reada(mcmcard,'DELE',dele,20.0d0)
5017 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
5018 call readi(mcmcard,'IREF',iref,0)
5019 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
5020 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
5021 call readi(mcmcard,'NCONF_IN',nconf_in,0)
5022 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
5023 write (iout,*) "NCONF_IN",nconf_in
5025 end subroutine csaread
5026 !-----------------------------------------------------------------------------
5030 use control_data, only: MaxMoveType
5033 ! implicit real*8 (a-h,o-z)
5034 ! include 'DIMENSIONS'
5035 ! include 'COMMON.MCM'
5036 ! include 'COMMON.MCE'
5037 ! include 'COMMON.IOUNITS'
5038 ! character(len=80) :: ucase
5039 character(len=320) :: mcmcard
5042 ! real(kind=8) :: var,ene
5044 call card_concat(mcmcard,.true.)
5045 call readi(mcmcard,'MAXACC',maxacc,100)
5046 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
5047 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
5048 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
5049 call readi(mcmcard,'MAXREPM',maxrepm,200)
5050 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
5051 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
5052 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
5053 call reada(mcmcard,'E_UP',e_up,5.0D0)
5054 call reada(mcmcard,'DELTE',delte,0.1D0)
5055 call readi(mcmcard,'NSWEEP',nsweep,5)
5056 call readi(mcmcard,'NSTEPH',nsteph,0)
5057 call readi(mcmcard,'NSTEPC',nstepc,0)
5058 call reada(mcmcard,'TMIN',tmin,298.0D0)
5059 call reada(mcmcard,'TMAX',tmax,298.0D0)
5060 call readi(mcmcard,'NWINDOW',nwindow,0)
5061 call readi(mcmcard,'PRINT_MC',print_mc,0)
5062 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
5063 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
5064 ent_read=(index(mcmcard,'ENT_READ').gt.0)
5065 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
5066 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
5067 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
5068 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
5069 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
5070 if (nwindow.gt.0) then
5071 allocate(winstart(nwindow)) !!el (maxres)
5072 allocate(winend(nwindow)) !!el
5073 allocate(winlen(nwindow)) !!el
5074 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
5076 winlen(i)=winend(i)-winstart(i)+1
5079 if (tmax.lt.tmin) tmax=tmin
5080 if (tmax.eq.tmin) then
5084 if (nstepc.gt.0 .and. nsteph.gt.0) then
5085 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
5086 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
5088 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
5089 ! Probabilities of different move types
5090 sumpro_type(0)=0.0D0
5091 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
5092 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
5093 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
5094 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
5095 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
5096 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
5097 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
5099 print *,'i',i,' sumprotype',sumpro_type(i)
5100 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
5101 print *,'i',i,' sumprotype',sumpro_type(i)
5104 end subroutine mcmread
5105 !-----------------------------------------------------------------------------
5106 subroutine read_minim
5109 ! implicit real*8 (a-h,o-z)
5110 ! include 'DIMENSIONS'
5111 ! include 'COMMON.MINIM'
5112 ! include 'COMMON.IOUNITS'
5113 ! character(len=80) :: ucase
5114 character(len=320) :: minimcard
5116 ! integer :: ntf,ik,iw_pdb
5117 ! real(kind=8) :: var,ene
5119 call card_concat(minimcard,.true.)
5120 call readi(minimcard,'MAXMIN',maxmin,2000)
5121 call readi(minimcard,'MAXFUN',maxfun,5000)
5122 call readi(minimcard,'MINMIN',minmin,maxmin)
5123 call readi(minimcard,'MINFUN',minfun,maxmin)
5124 call reada(minimcard,'TOLF',tolf,1.0D-2)
5125 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
5126 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
5127 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
5128 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
5129 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
5130 'Options in energy minimization:'
5131 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
5132 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
5133 'MinMin:',MinMin,' MinFun:',MinFun,&
5134 ' TolF:',TolF,' RTolF:',RTolF
5136 end subroutine read_minim
5137 !-----------------------------------------------------------------------------
5138 subroutine openunits
5140 use MD_data, only: usampl
5143 use control_data, only:out1file
5144 use control, only: getenv_loc
5145 ! implicit real*8 (a-h,o-z)
5146 ! include 'DIMENSIONS'
5149 character(len=16) :: form,nodename
5150 integer :: nodelen,ierror,npos
5152 ! include 'COMMON.SETUP'
5153 ! include 'COMMON.IOUNITS'
5154 ! include 'COMMON.MD'
5155 ! include 'COMMON.CONTROL'
5156 integer :: lenpre,lenpot,lentmp !,ilen
5158 character(len=3) :: out1file_text !,ucase
5159 character(len=3) :: ll
5162 ! integer :: ntf,ik,iw_pdb
5163 ! real(kind=8) :: var,ene
5165 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
5166 call getenv_loc("PREFIX",prefix)
5168 call getenv_loc("POT",pot)
5169 call getenv_loc("DIRTMP",tmpdir)
5170 call getenv_loc("CURDIR",curdir)
5171 call getenv_loc("OUT1FILE",out1file_text)
5172 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
5173 out1file_text=ucase(out1file_text)
5174 if (out1file_text(1:1).eq."Y") then
5177 out1file=fg_rank.gt.0
5182 if (lentmp.gt.0) then
5183 write (*,'(80(1h!))')
5184 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
5185 write (*,'(80(1h!))')
5186 write (*,*)"All output files will be on node /tmp directory."
5188 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
5189 if (me.eq.king) then
5190 write (*,*) "The master node is ",nodename
5191 else if (fg_rank.eq.0) then
5192 write (*,*) "I am the CG slave node ",nodename
5194 write (*,*) "I am the FG slave node ",nodename
5197 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
5198 lenpre = lentmp+lenpre+1
5200 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
5201 ! Get the names and open the input files
5202 #if defined(WINIFL) || defined(WINPGI)
5203 open(1,file=pref_orig(:ilen(pref_orig))// &
5204 '.inp',status='old',readonly,shared)
5205 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5206 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5207 ! Get parameter filenames and open the parameter files.
5208 call getenv_loc('BONDPAR',bondname)
5209 open (ibond,file=bondname,status='old',readonly,shared)
5210 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5211 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
5212 call getenv_loc('THETPAR',thetname)
5213 open (ithep,file=thetname,status='old',readonly,shared)
5214 call getenv_loc('ROTPAR',rotname)
5215 open (irotam,file=rotname,status='old',readonly,shared)
5216 call getenv_loc('TORPAR',torname)
5217 open (itorp,file=torname,status='old',readonly,shared)
5218 call getenv_loc('TORDPAR',tordname)
5219 open (itordp,file=tordname,status='old',readonly,shared)
5220 call getenv_loc('FOURIER',fouriername)
5221 open (ifourier,file=fouriername,status='old',readonly,shared)
5222 call getenv_loc('ELEPAR',elename)
5223 open (ielep,file=elename,status='old',readonly,shared)
5224 call getenv_loc('SIDEPAR',sidename)
5225 open (isidep,file=sidename,status='old',readonly,shared)
5227 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5228 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
5229 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5230 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
5231 call getenv_loc('TORPAR_NUCL',torname_nucl)
5232 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
5233 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5234 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
5235 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5236 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
5239 #elif (defined CRAY) || (defined AIX)
5240 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5242 ! print *,"Processor",myrank," opened file 1"
5243 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5244 ! print *,"Processor",myrank," opened file 9"
5245 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5246 ! Get parameter filenames and open the parameter files.
5247 call getenv_loc('BONDPAR',bondname)
5248 open (ibond,file=bondname,status='old',action='read')
5249 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5250 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5252 ! print *,"Processor",myrank," opened file IBOND"
5253 call getenv_loc('THETPAR',thetname)
5254 open (ithep,file=thetname,status='old',action='read')
5255 ! print *,"Processor",myrank," opened file ITHEP"
5256 call getenv_loc('ROTPAR',rotname)
5257 open (irotam,file=rotname,status='old',action='read')
5258 ! print *,"Processor",myrank," opened file IROTAM"
5259 call getenv_loc('TORPAR',torname)
5260 open (itorp,file=torname,status='old',action='read')
5261 ! print *,"Processor",myrank," opened file ITORP"
5262 call getenv_loc('TORDPAR',tordname)
5263 open (itordp,file=tordname,status='old',action='read')
5264 ! print *,"Processor",myrank," opened file ITORDP"
5265 call getenv_loc('SCCORPAR',sccorname)
5266 open (isccor,file=sccorname,status='old',action='read')
5267 ! print *,"Processor",myrank," opened file ISCCOR"
5268 call getenv_loc('FOURIER',fouriername)
5269 open (ifourier,file=fouriername,status='old',action='read')
5270 ! print *,"Processor",myrank," opened file IFOURIER"
5271 call getenv_loc('ELEPAR',elename)
5272 open (ielep,file=elename,status='old',action='read')
5273 ! print *,"Processor",myrank," opened file IELEP"
5274 call getenv_loc('SIDEPAR',sidename)
5275 open (isidep,file=sidename,status='old',action='read')
5277 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5278 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5279 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5280 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5281 call getenv_loc('TORPAR_NUCL',torname_nucl)
5282 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5283 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5284 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5285 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5286 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5288 call getenv_loc('LIPTRANPAR',liptranname)
5289 open (iliptranpar,file=liptranname,status='old',action='read')
5290 call getenv_loc('TUBEPAR',tubename)
5291 open (itube,file=tubename,status='old',action='read')
5292 call getenv_loc('IONPAR',ionname)
5293 open (iion,file=ionname,status='old',action='read')
5295 ! print *,"Processor",myrank," opened file ISIDEP"
5296 ! print *,"Processor",myrank," opened parameter files"
5298 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
5299 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5300 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5301 ! Get parameter filenames and open the parameter files.
5302 call getenv_loc('BONDPAR',bondname)
5303 open (ibond,file=bondname,status='old')
5304 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5305 open (ibond_nucl,file=bondname_nucl,status='old')
5307 call getenv_loc('THETPAR',thetname)
5308 open (ithep,file=thetname,status='old')
5309 call getenv_loc('ROTPAR',rotname)
5310 open (irotam,file=rotname,status='old')
5311 call getenv_loc('TORPAR',torname)
5312 open (itorp,file=torname,status='old')
5313 call getenv_loc('TORDPAR',tordname)
5314 open (itordp,file=tordname,status='old')
5315 call getenv_loc('SCCORPAR',sccorname)
5316 open (isccor,file=sccorname,status='old')
5317 call getenv_loc('FOURIER',fouriername)
5318 open (ifourier,file=fouriername,status='old')
5319 call getenv_loc('ELEPAR',elename)
5320 open (ielep,file=elename,status='old')
5321 call getenv_loc('SIDEPAR',sidename)
5322 open (isidep,file=sidename,status='old')
5324 open (ithep_nucl,file=thetname_nucl,status='old')
5325 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5326 open (irotam_nucl,file=rotname_nucl,status='old')
5327 call getenv_loc('TORPAR_NUCL',torname_nucl)
5328 open (itorp_nucl,file=torname_nucl,status='old')
5329 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5330 open (itordp_nucl,file=tordname_nucl,status='old')
5331 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5332 open (isidep_nucl,file=sidename_nucl,status='old')
5334 call getenv_loc('LIPTRANPAR',liptranname)
5335 open (iliptranpar,file=liptranname,status='old')
5336 call getenv_loc('TUBEPAR',tubename)
5337 open (itube,file=tubename,status='old')
5338 call getenv_loc('IONPAR',ionname)
5339 open (iion,file=ionname,status='old')
5341 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5343 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5344 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5345 ! Get parameter filenames and open the parameter files.
5346 call getenv_loc('BONDPAR',bondname)
5347 open (ibond,file=bondname,status='old',action='read')
5348 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5349 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5350 call getenv_loc('THETPAR',thetname)
5351 open (ithep,file=thetname,status='old',action='read')
5352 call getenv_loc('ROTPAR',rotname)
5353 open (irotam,file=rotname,status='old',action='read')
5354 call getenv_loc('TORPAR',torname)
5355 open (itorp,file=torname,status='old',action='read')
5356 call getenv_loc('TORDPAR',tordname)
5357 open (itordp,file=tordname,status='old',action='read')
5358 call getenv_loc('SCCORPAR',sccorname)
5359 open (isccor,file=sccorname,status='old',action='read')
5361 call getenv_loc('THETPARPDB',thetname_pdb)
5362 print *,"thetname_pdb ",thetname_pdb
5363 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
5364 print *,ithep_pdb," opened"
5366 call getenv_loc('FOURIER',fouriername)
5367 open (ifourier,file=fouriername,status='old',readonly)
5368 call getenv_loc('ELEPAR',elename)
5369 open (ielep,file=elename,status='old',readonly)
5370 call getenv_loc('SIDEPAR',sidename)
5371 open (isidep,file=sidename,status='old',readonly)
5373 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5374 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5375 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5376 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5377 call getenv_loc('TORPAR_NUCL',torname_nucl)
5378 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5379 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5380 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5381 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5382 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5383 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
5384 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
5385 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
5386 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
5387 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
5388 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
5389 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
5390 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
5393 call getenv_loc('LIPTRANPAR',liptranname)
5394 open (iliptranpar,file=liptranname,status='old',action='read')
5395 call getenv_loc('TUBEPAR',tubename)
5396 open (itube,file=tubename,status='old',action='read')
5397 call getenv_loc('IONPAR',ionname)
5398 open (iion,file=ionname,status='old',action='read')
5401 call getenv_loc('ROTPARPDB',rotname_pdb)
5402 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
5405 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
5406 #if defined(WINIFL) || defined(WINPGI)
5407 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
5408 #elif (defined CRAY) || (defined AIX)
5409 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5411 open (iscpp_nucl,file=scpname_nucl,status='old')
5413 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5418 ! 8/9/01 In the newest version SCp interaction constants are read from a file
5419 ! Use -DOLDSCP to use hard-coded constants instead.
5421 call getenv_loc('SCPPAR',scpname)
5422 #if defined(WINIFL) || defined(WINPGI)
5423 open (iscpp,file=scpname,status='old',readonly,shared)
5424 #elif (defined CRAY) || (defined AIX)
5425 open (iscpp,file=scpname,status='old',action='read')
5427 open (iscpp,file=scpname,status='old')
5429 open (iscpp,file=scpname,status='old',action='read')
5432 call getenv_loc('PATTERN',patname)
5433 #if defined(WINIFL) || defined(WINPGI)
5434 open (icbase,file=patname,status='old',readonly,shared)
5435 #elif (defined CRAY) || (defined AIX)
5436 open (icbase,file=patname,status='old',action='read')
5438 open (icbase,file=patname,status='old')
5440 open (icbase,file=patname,status='old',action='read')
5443 ! Open output file only for CG processes
5444 ! print *,"Processor",myrank," fg_rank",fg_rank
5445 if (fg_rank.eq.0) then
5447 if (nodes.eq.1) then
5450 npos = dlog10(dfloat(nodes-1))+1
5452 if (npos.lt.3) npos=3
5453 write (liczba,'(i1)') npos
5454 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
5456 write (liczba,form) me
5457 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
5458 liczba(:ilen(liczba))
5459 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5461 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5463 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
5464 liczba(:ilen(liczba))//'.mol2'
5465 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5466 liczba(:ilen(liczba))//'.stat'
5468 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
5469 //liczba(:ilen(liczba))//'.stat')
5470 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
5473 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5474 liczba(:ilen(liczba))//'.const'
5479 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
5480 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
5481 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
5482 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
5483 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
5485 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
5487 rest2name=prefix(:ilen(prefix))//'.rst'
5489 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
5492 #if defined(AIX) || defined(PGI)
5493 if (me.eq.king .or. .not. out1file) &
5494 open(iout,file=outname,status='unknown')
5496 if (fg_rank.gt.0) then
5497 write (liczba,'(i3.3)') myrank/nfgtasks
5498 write (ll,'(bz,i3.3)') fg_rank
5499 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5504 open(igeom,file=intname,status='unknown',position='append')
5505 open(ipdb,file=pdbname,status='unknown')
5506 open(imol2,file=mol2name,status='unknown')
5507 open(istat,file=statname,status='unknown',position='append')
5509 !1out open(iout,file=outname,status='unknown')
5512 if (me.eq.king .or. .not.out1file) &
5513 open(iout,file=outname,status='unknown')
5515 if (fg_rank.gt.0) then
5516 write (liczba,'(i3.3)') myrank/nfgtasks
5517 write (ll,'(bz,i3.3)') fg_rank
5518 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5523 open(igeom,file=intname,status='unknown',access='append')
5524 open(ipdb,file=pdbname,status='unknown')
5525 open(imol2,file=mol2name,status='unknown')
5526 open(istat,file=statname,status='unknown',access='append')
5528 !1out open(iout,file=outname,status='unknown')
5531 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
5532 csa_seed=prefix(:lenpre)//'.CSA.seed'
5533 csa_history=prefix(:lenpre)//'.CSA.history'
5534 csa_bank=prefix(:lenpre)//'.CSA.bank'
5535 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
5536 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
5537 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
5538 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
5539 csa_int=prefix(:lenpre)//'.int'
5540 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
5541 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
5542 csa_in=prefix(:lenpre)//'.CSA.in'
5543 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
5546 write (iout,'(80(1h-))')
5547 write (iout,'(30x,a)') "FILE ASSIGNMENT"
5548 write (iout,'(80(1h-))')
5549 write (iout,*) "Input file : ",&
5550 pref_orig(:ilen(pref_orig))//'.inp'
5551 write (iout,*) "Output file : ",&
5552 outname(:ilen(outname))
5554 write (iout,*) "Sidechain potential file : ",&
5555 sidename(:ilen(sidename))
5557 write (iout,*) "SCp potential file : ",&
5558 scpname(:ilen(scpname))
5560 write (iout,*) "Electrostatic potential file : ",&
5561 elename(:ilen(elename))
5562 write (iout,*) "Cumulant coefficient file : ",&
5563 fouriername(:ilen(fouriername))
5564 write (iout,*) "Torsional parameter file : ",&
5565 torname(:ilen(torname))
5566 write (iout,*) "Double torsional parameter file : ",&
5567 tordname(:ilen(tordname))
5568 write (iout,*) "SCCOR parameter file : ",&
5569 sccorname(:ilen(sccorname))
5570 write (iout,*) "Bond & inertia constant file : ",&
5571 bondname(:ilen(bondname))
5572 write (iout,*) "Bending parameter file : ",&
5573 thetname(:ilen(thetname))
5574 write (iout,*) "Rotamer parameter file : ",&
5575 rotname(:ilen(rotname))
5578 write (iout,*) "Thetpdb parameter file : ",&
5579 thetname_pdb(:ilen(thetname_pdb))
5582 write (iout,*) "Threading database : ",&
5583 patname(:ilen(patname))
5585 write (iout,*)" DIRTMP : ",&
5587 write (iout,'(80(1h-))')
5590 end subroutine openunits
5591 !-----------------------------------------------------------------------------
5594 use geometry_data, only: nres,dc
5596 ! implicit real*8 (a-h,o-z)
5597 ! include 'DIMENSIONS'
5598 ! include 'COMMON.CHAIN'
5599 ! include 'COMMON.IOUNITS'
5600 ! include 'COMMON.MD'
5603 ! real(kind=8) :: var,ene
5605 open(irest2,file=rest2name,status='unknown')
5606 read(irest2,*) totT,EK,potE,totE,t_bath
5609 ! AL 4/17/17: Now reading d_t(0,:) too
5611 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
5614 ! AL 4/17/17: Now reading d_c(0,:) too
5616 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
5619 read (irest2,*) iset
5623 end subroutine readrst
5624 !-----------------------------------------------------------------------------
5625 subroutine read_fragments
5629 use control_data, only:out1file
5632 ! implicit real*8 (a-h,o-z)
5633 ! include 'DIMENSIONS'
5637 ! include 'COMMON.SETUP'
5638 ! include 'COMMON.CHAIN'
5639 ! include 'COMMON.IOUNITS'
5640 ! include 'COMMON.MD'
5641 ! include 'COMMON.CONTROL'
5644 ! real(kind=8) :: var,ene
5646 read(inp,*) nset,nfrag,npair,nfrag_back
5648 !el from module energy
5649 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
5650 if(.not.allocated(wfrag_back)) then
5651 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5652 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5654 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
5655 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
5657 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
5658 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
5661 if(me.eq.king.or..not.out1file) &
5662 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
5663 " nfrag_back",nfrag_back
5665 read(inp,*) mset(iset)
5667 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
5669 if(me.eq.king.or..not.out1file) &
5670 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
5671 ifrag(2,i,iset), qinfrag(i,iset)
5674 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
5676 if(me.eq.king.or..not.out1file) &
5677 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
5678 ipair(2,i,iset), qinpair(i,iset)
5681 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5682 wfrag_back(3,i,iset),&
5683 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5684 if(me.eq.king.or..not.out1file) &
5685 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5686 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5690 end subroutine read_fragments
5691 !-----------------------------------------------------------------------------
5693 !-----------------------------------------------------------------------------
5697 ! implicit real*8 (a-h,o-z)
5698 ! include 'DIMENSIONS'
5699 ! include 'COMMON.CSA'
5700 ! include 'COMMON.BANK'
5701 ! include 'COMMON.IOUNITS'
5703 ! integer :: ntf,ik,iw_pdb
5704 ! real(kind=8) :: var,ene
5706 open(icsa_in,file=csa_in,status="old",err=100)
5707 read(icsa_in,*) nconf
5708 read(icsa_in,*) jstart,jend
5709 read(icsa_in,*) nstmax
5710 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5711 read(icsa_in,*) nran0,nran1,irr
5712 read(icsa_in,*) nseed
5713 read(icsa_in,*) ntotal,cut1,cut2
5714 read(icsa_in,*) estop
5715 read(icsa_in,*) icmax,irestart
5716 read(icsa_in,*) ntbankm,dele,difcut
5717 read(icsa_in,*) iref,rmscut,pnccut
5718 read(icsa_in,*) ndiff
5725 end subroutine csa_read
5726 !-----------------------------------------------------------------------------
5727 subroutine initial_write
5730 ! implicit real*8 (a-h,o-z)
5731 ! include 'DIMENSIONS'
5732 ! include 'COMMON.CSA'
5733 ! include 'COMMON.BANK'
5734 ! include 'COMMON.IOUNITS'
5736 ! integer :: ntf,ik,iw_pdb
5737 ! real(kind=8) :: var,ene
5739 open(icsa_seed,file=csa_seed,status="unknown")
5740 write(icsa_seed,*) "seed"
5742 #if defined(AIX) || defined(PGI)
5743 open(icsa_history,file=csa_history,status="unknown",&
5746 open(icsa_history,file=csa_history,status="unknown",&
5749 write(icsa_history,*) nconf
5750 write(icsa_history,*) jstart,jend
5751 write(icsa_history,*) nstmax
5752 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5753 write(icsa_history,*) nran0,nran1,irr
5754 write(icsa_history,*) nseed
5755 write(icsa_history,*) ntotal,cut1,cut2
5756 write(icsa_history,*) estop
5757 write(icsa_history,*) icmax,irestart
5758 write(icsa_history,*) ntbankm,dele,difcut
5759 write(icsa_history,*) iref,rmscut,pnccut
5760 write(icsa_history,*) ndiff
5762 write(icsa_history,*)
5765 open(icsa_bank1,file=csa_bank1,status="unknown")
5766 write(icsa_bank1,*) 0
5770 end subroutine initial_write
5771 !-----------------------------------------------------------------------------
5772 subroutine restart_write
5775 ! implicit real*8 (a-h,o-z)
5776 ! include 'DIMENSIONS'
5777 ! include 'COMMON.IOUNITS'
5778 ! include 'COMMON.CSA'
5779 ! include 'COMMON.BANK'
5781 ! integer :: ntf,ik,iw_pdb
5782 ! real(kind=8) :: var,ene
5784 #if defined(AIX) || defined(PGI)
5785 open(icsa_history,file=csa_history,position="append")
5787 open(icsa_history,file=csa_history,access="append")
5789 write(icsa_history,*)
5790 write(icsa_history,*) "This is restart"
5791 write(icsa_history,*)
5792 write(icsa_history,*) nconf
5793 write(icsa_history,*) jstart,jend
5794 write(icsa_history,*) nstmax
5795 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5796 write(icsa_history,*) nran0,nran1,irr
5797 write(icsa_history,*) nseed
5798 write(icsa_history,*) ntotal,cut1,cut2
5799 write(icsa_history,*) estop
5800 write(icsa_history,*) icmax,irestart
5801 write(icsa_history,*) ntbankm,dele,difcut
5802 write(icsa_history,*) iref,rmscut,pnccut
5803 write(icsa_history,*) ndiff
5804 write(icsa_history,*)
5805 write(icsa_history,*) "irestart is: ", irestart
5807 write(icsa_history,*)
5811 end subroutine restart_write
5812 !-----------------------------------------------------------------------------
5814 !-----------------------------------------------------------------------------
5815 subroutine write_pdb(npdb,titelloc,ee)
5817 ! implicit real*8 (a-h,o-z)
5818 ! include 'DIMENSIONS'
5819 ! include 'COMMON.IOUNITS'
5820 character(len=50) :: titelloc1
5821 character*(*) :: titelloc
5822 character(len=3) :: zahl
5823 character(len=5) :: liczba5
5825 integer :: npdb !,ilen
5829 ! real(kind=8) :: var,ene
5833 if (npdb.lt.1000) then
5834 call numstr(npdb,zahl)
5835 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5837 if (npdb.lt.10000) then
5838 write(liczba5,'(i1,i4)') 0,npdb
5840 write(liczba5,'(i5)') npdb
5842 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5844 call pdbout(ee,titelloc1,ipdb)
5847 end subroutine write_pdb
5848 !-----------------------------------------------------------------------------
5850 !-----------------------------------------------------------------------------
5851 subroutine write_thread_summary
5852 ! Thread the sequence through a database of known structures
5853 use control_data, only: refstr
5855 use energy_data, only: n_ene_comp
5857 ! implicit real*8 (a-h,o-z)
5858 ! include 'DIMENSIONS'
5860 use MPI_data !include 'COMMON.INFO'
5863 ! include 'COMMON.CONTROL'
5864 ! include 'COMMON.CHAIN'
5865 ! include 'COMMON.DBASE'
5866 ! include 'COMMON.INTERACT'
5867 ! include 'COMMON.VAR'
5868 ! include 'COMMON.THREAD'
5869 ! include 'COMMON.FFIELD'
5870 ! include 'COMMON.SBRIDGE'
5871 ! include 'COMMON.HEADER'
5872 ! include 'COMMON.NAMES'
5873 ! include 'COMMON.IOUNITS'
5874 ! include 'COMMON.TIME1'
5876 integer,dimension(maxthread) :: ip
5877 real(kind=8),dimension(0:n_ene) :: energia
5879 integer :: i,j,ii,jj,ipj,ik,kk,ist
5880 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5882 write (iout,'(30x,a/)') &
5883 ' *********** Summary threading statistics ************'
5884 write (iout,'(a)') 'Initial energies:'
5885 write (iout,'(a4,2x,a12,14a14,3a8)') &
5886 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5887 'RMSnat','NatCONT','NNCONT','RMS'
5888 ! Energy sort patterns
5893 enet=ener(n_ene-1,ip(i))
5896 if (ener(n_ene-1,ip(j)).lt.enet) then
5898 enet=ener(n_ene-1,ip(j))
5910 ist=nres_base(2,ii)+ipatt(2,i)
5912 energia(i)=ener0(kk,i)
5914 etot=ener0(n_ene_comp+1,i)
5915 rmsnat=ener0(n_ene_comp+2,i)
5916 rms=ener0(n_ene_comp+3,i)
5917 frac=ener0(n_ene_comp+4,i)
5918 frac_nn=ener0(n_ene_comp+5,i)
5921 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5922 i,str_nam(ii),ist+1,&
5923 (energia(print_order(kk)),kk=1,nprint_ene),&
5924 etot,rmsnat,frac,frac_nn,rms
5926 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5927 i,str_nam(ii),ist+1,&
5928 (energia(print_order(kk)),kk=1,nprint_ene),etot
5931 write (iout,'(//a)') 'Final energies:'
5932 write (iout,'(a4,2x,a12,17a14,3a8)') &
5933 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5934 'RMSnat','NatCONT','NNCONT','RMS'
5938 ist=nres_base(2,ii)+ipatt(2,i)
5940 energia(kk)=ener(kk,ik)
5942 etot=ener(n_ene_comp+1,i)
5943 rmsnat=ener(n_ene_comp+2,i)
5944 rms=ener(n_ene_comp+3,i)
5945 frac=ener(n_ene_comp+4,i)
5946 frac_nn=ener(n_ene_comp+5,i)
5947 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5948 i,str_nam(ii),ist+1,&
5949 (energia(print_order(kk)),kk=1,nprint_ene),&
5950 etot,rmsnat,frac,frac_nn,rms
5952 write (iout,'(/a/)') 'IEXAM array:'
5953 write (iout,'(i5)') nexcl
5955 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5957 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5958 'Max. time for threading step ',max_time_for_thread,&
5959 'Average time for threading step: ',ave_time_for_thread
5961 end subroutine write_thread_summary
5962 !-----------------------------------------------------------------------------
5963 subroutine write_stat_thread(ithread,ipattern,ist)
5965 use energy_data, only: n_ene_comp
5967 ! implicit real*8 (a-h,o-z)
5968 ! include "DIMENSIONS"
5969 ! include "COMMON.CONTROL"
5970 ! include "COMMON.IOUNITS"
5971 ! include "COMMON.THREAD"
5972 ! include "COMMON.FFIELD"
5973 ! include "COMMON.DBASE"
5974 ! include "COMMON.NAMES"
5975 real(kind=8),dimension(0:n_ene) :: energia
5977 integer :: ithread,ipattern,ist,i
5978 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5980 #if defined(AIX) || defined(PGI)
5981 open(istat,file=statname,position='append')
5983 open(istat,file=statname,access='append')
5986 energia(i)=ener(i,ithread)
5988 etot=ener(n_ene_comp+1,ithread)
5989 rmsnat=ener(n_ene_comp+2,ithread)
5990 rms=ener(n_ene_comp+3,ithread)
5991 frac=ener(n_ene_comp+4,ithread)
5992 frac_nn=ener(n_ene_comp+5,ithread)
5993 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5994 ithread,str_nam(ipattern),ist+1,&
5995 (energia(print_order(i)),i=1,nprint_ene),&
5996 etot,rmsnat,frac,frac_nn,rms
5999 end subroutine write_stat_thread
6000 !-----------------------------------------------------------------------------
6002 !-----------------------------------------------------------------------------
6003 end module io_config