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
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
743 integer :: maxinter,junk,kk,ii,ncatprotparm
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
750 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
751 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
752 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 ! For printing parameters after they are read set the following in the UNRES
757 ! setenv PRINT_PARM YES
759 ! To print parameters in LaTeX format rather than as ASCII tables:
763 call getenv_loc("PRINT_PARM",lancuch)
764 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
765 call getenv_loc("LATEX",lancuch)
766 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
768 dwa16=2.0d0**(1.0d0/6.0d0)
770 ! Assign virtual-bond length
773 vblinv2=vblinv*vblinv
775 ! Read the virtual-bond parameters, masses, and moments of inertia
776 ! and Stokes' radii of the peptide group and side chains
778 allocate(dsc(ntyp1)) !(ntyp1)
779 allocate(dsc_inv(ntyp1)) !(ntyp1)
780 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
781 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
782 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
783 allocate(nbondterm(ntyp)) !(ntyp)
784 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
785 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
786 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(long_r_sidechain(ntyp))
788 allocate(short_r_sidechain(ntyp))
793 allocate(msc(ntyp+1)) !(ntyp+1)
794 allocate(isc(ntyp+1)) !(ntyp+1)
795 allocate(restok(ntyp+1)) !(ntyp+1)
797 read (ibond,*) vbldp0,akp,mp,ip,pstok
800 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
801 dsc(i) = vbldsc0(1,i)
805 dsc_inv(i)=1.0D0/dsc(i)
814 allocate(msc(ntyp+1,5)) !(ntyp+1)
815 allocate(isc(ntyp+1,5)) !(ntyp+1)
816 allocate(restok(ntyp+1,5)) !(ntyp+1)
818 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
820 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
821 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
822 dsc(i) = vbldsc0(1,i)
826 dsc_inv(i)=1.0D0/dsc(i)
830 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
833 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
834 ! dsc(i) = vbldsc0_nucl(1,i)
838 ! dsc_inv(i)=1.0D0/dsc(i)
841 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
842 ! do i=1,ntyp_molec(2)
843 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
844 ! aksc_nucl(j,i),abond0_nucl(j,i),&
845 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
846 ! dsc(i) = vbldsc0(1,i)
850 ! dsc_inv(i)=1.0D0/dsc(i)
855 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
856 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
858 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
860 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
861 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
863 write (iout,'(13x,3f10.5)') &
864 vbldsc0(j,i),aksc(j,i),abond0(j,i)
869 read(iion,*) msc(i,5),restok(i,5)
870 print *,msc(i,5),restok(i,5)
874 read (iion,*) ncatprotparm
875 allocate(catprm(ncatprotparm,4))
877 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
880 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
881 !----------------------------------------------------
882 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
883 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
884 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
885 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
886 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
887 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
897 allocate(liptranene(ntyp))
898 !C reading lipid parameters
899 write (iout,*) "iliptranpar",iliptranpar
901 read(iliptranpar,*) pepliptran
904 read(iliptranpar,*) liptranene(i)
905 print *,liptranene(i)
911 ! Read the parameters of the probability distribution/energy expression
912 ! of the virtual-bond valence angles theta
915 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
916 (bthet(j,i,1,1),j=1,2)
917 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
918 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
919 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
923 athet(1,i,1,-1)=athet(1,i,1,1)
924 athet(2,i,1,-1)=athet(2,i,1,1)
925 bthet(1,i,1,-1)=-bthet(1,i,1,1)
926 bthet(2,i,1,-1)=-bthet(2,i,1,1)
927 athet(1,i,-1,1)=-athet(1,i,1,1)
928 athet(2,i,-1,1)=-athet(2,i,1,1)
929 bthet(1,i,-1,1)=bthet(1,i,1,1)
930 bthet(2,i,-1,1)=bthet(2,i,1,1)
934 athet(1,i,-1,-1)=athet(1,-i,1,1)
935 athet(2,i,-1,-1)=-athet(2,-i,1,1)
936 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
937 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
938 athet(1,i,-1,1)=athet(1,-i,1,1)
939 athet(2,i,-1,1)=-athet(2,-i,1,1)
940 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
941 bthet(2,i,-1,1)=bthet(2,-i,1,1)
942 athet(1,i,1,-1)=-athet(1,-i,1,1)
943 athet(2,i,1,-1)=athet(2,-i,1,1)
944 bthet(1,i,1,-1)=bthet(1,-i,1,1)
945 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
950 polthet(j,i)=polthet(j,-i)
953 gthet(j,i)=gthet(j,-i)
961 'Parameters of the virtual-bond valence angles:'
962 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
963 ' ATHETA0 ',' A1 ',' A2 ',&
966 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
967 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
969 write (iout,'(/a/9x,5a/79(1h-))') &
970 'Parameters of the expression for sigma(theta_c):',&
971 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
972 ' ALPH3 ',' SIGMA0C '
974 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
975 (polthet(j,i),j=0,3),sigc0(i)
977 write (iout,'(/a/9x,5a/79(1h-))') &
978 'Parameters of the second gaussian:',&
979 ' THETA0 ',' SIGMA0 ',' G1 ',&
982 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
983 sig0(i),(gthet(j,i),j=1,3)
987 'Parameters of the virtual-bond valence angles:'
988 write (iout,'(/a/9x,5a/79(1h-))') &
989 'Coefficients of expansion',&
990 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
991 ' b1*10^1 ',' b2*10^1 '
993 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
994 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
995 (10*bthet(j,i,1,1),j=1,2)
997 write (iout,'(/a/9x,5a/79(1h-))') &
998 'Parameters of the expression for sigma(theta_c):',&
999 ' alpha0 ',' alph1 ',' alph2 ',&
1000 ' alhp3 ',' sigma0c '
1002 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1003 (polthet(j,i),j=0,3),sigc0(i)
1005 write (iout,'(/a/9x,5a/79(1h-))') &
1006 'Parameters of the second gaussian:',&
1007 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1010 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1011 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1017 ! Read the parameters of Utheta determined from ab initio surfaces
1018 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1020 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1021 ntheterm3,nsingle,ndouble
1022 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1024 !----------------------------------------------------
1025 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1026 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1027 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1028 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1029 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1030 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1031 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1032 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1033 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1034 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1035 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1036 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1037 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1038 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1039 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1040 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1041 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1042 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1043 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1044 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1045 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1046 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1047 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1048 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1050 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1052 ithetyp(i)=-ithetyp(-i)
1055 aa0thet(:,:,:,:)=0.0d0
1056 aathet(:,:,:,:,:)=0.0d0
1057 bbthet(:,:,:,:,:,:)=0.0d0
1058 ccthet(:,:,:,:,:,:)=0.0d0
1059 ddthet(:,:,:,:,:,:)=0.0d0
1060 eethet(:,:,:,:,:,:)=0.0d0
1061 ffthet(:,:,:,:,:,:,:)=0.0d0
1062 ggthet(:,:,:,:,:,:,:)=0.0d0
1064 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1066 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1067 ! VAR:1=non-glicyne non-proline 2=proline
1068 ! VAR:negative values for D-aminoacid
1070 do j=-nthetyp,nthetyp
1071 do k=-nthetyp,nthetyp
1072 read (ithep,'(6a)',end=111,err=111) res1
1073 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1074 ! VAR: aa0thet is variable describing the average value of Foureir
1075 ! VAR: expansion series
1076 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1077 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1078 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1079 read (ithep,*,end=111,err=111) &
1080 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1081 read (ithep,*,end=111,err=111) &
1082 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1083 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1084 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1085 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1087 read (ithep,*,end=111,err=111) &
1088 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1089 ffthet(lll,llll,ll,i,j,k,iblock),&
1090 ggthet(llll,lll,ll,i,j,k,iblock),&
1091 ggthet(lll,llll,ll,i,j,k,iblock),&
1092 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1097 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1098 ! coefficients of theta-and-gamma-dependent terms are zero.
1099 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1100 ! RECOMENTDED AFTER VERSION 3.3)
1104 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1105 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1107 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1108 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1111 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1113 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1116 ! AND COMMENT THE LOOPS BELOW
1120 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1121 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1123 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1124 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1127 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1129 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1134 ! Substitution for D aminoacids from symmetry.
1137 do j=-nthetyp,nthetyp
1138 do k=-nthetyp,nthetyp
1139 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1141 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1145 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1146 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1147 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1148 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1154 ffthet(llll,lll,ll,i,j,k,iblock)= &
1155 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1156 ffthet(lll,llll,ll,i,j,k,iblock)= &
1157 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1158 ggthet(llll,lll,ll,i,j,k,iblock)= &
1159 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1160 ggthet(lll,llll,ll,i,j,k,iblock)= &
1161 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1170 ! Control printout of the coefficients of virtual-bond-angle potentials
1173 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1178 write (iout,'(//4a)') &
1179 'Type ',onelett(i),onelett(j),onelett(k)
1180 write (iout,'(//a,10x,a)') " l","a[l]"
1181 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1182 write (iout,'(i2,1pe15.5)') &
1183 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1185 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1186 "b",l,"c",l,"d",l,"e",l
1188 write (iout,'(i2,4(1pe15.5))') m,&
1189 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1190 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1194 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1195 "f+",l,"f-",l,"g+",l,"g-",l
1198 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1199 ffthet(n,m,l,i,j,k,iblock),&
1200 ffthet(m,n,l,i,j,k,iblock),&
1201 ggthet(n,m,l,i,j,k,iblock),&
1202 ggthet(m,n,l,i,j,k,iblock)
1212 write (2,*) "Start reading THETA_PDB",ithep_pdb
1214 ! write (2,*) 'i=',i
1215 read (ithep_pdb,*,err=111,end=111) &
1216 a0thet(i),(athet(j,i,1,1),j=1,2),&
1217 (bthet(j,i,1,1),j=1,2)
1218 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1219 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1220 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1221 sigc0(i)=sigc0(i)**2
1224 athet(1,i,1,-1)=athet(1,i,1,1)
1225 athet(2,i,1,-1)=athet(2,i,1,1)
1226 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1227 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1228 athet(1,i,-1,1)=-athet(1,i,1,1)
1229 athet(2,i,-1,1)=-athet(2,i,1,1)
1230 bthet(1,i,-1,1)=bthet(1,i,1,1)
1231 bthet(2,i,-1,1)=bthet(2,i,1,1)
1234 a0thet(i)=a0thet(-i)
1235 athet(1,i,-1,-1)=athet(1,-i,1,1)
1236 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1237 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1238 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1239 athet(1,i,-1,1)=athet(1,-i,1,1)
1240 athet(2,i,-1,1)=-athet(2,-i,1,1)
1241 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1242 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1243 athet(1,i,1,-1)=-athet(1,-i,1,1)
1244 athet(2,i,1,-1)=athet(2,-i,1,1)
1245 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1246 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1247 theta0(i)=theta0(-i)
1251 polthet(j,i)=polthet(j,-i)
1254 gthet(j,i)=gthet(j,-i)
1257 write (2,*) "End reading THETA_PDB"
1261 !--------------- Reading theta parameters for nucleic acid-------
1262 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1263 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1264 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1265 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1266 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1267 nthetyp_nucl+1,nthetyp_nucl+1))
1268 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1269 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1270 nthetyp_nucl+1,nthetyp_nucl+1))
1271 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1272 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1273 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1274 nthetyp_nucl+1,nthetyp_nucl+1))
1275 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1276 nthetyp_nucl+1,nthetyp_nucl+1))
1277 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1278 nthetyp_nucl+1,nthetyp_nucl+1))
1279 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1280 nthetyp_nucl+1,nthetyp_nucl+1))
1281 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1282 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1283 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1284 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1285 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1286 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1288 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1289 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1291 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1293 aa0thet_nucl(:,:,:)=0.0d0
1294 aathet_nucl(:,:,:,:)=0.0d0
1295 bbthet_nucl(:,:,:,:,:)=0.0d0
1296 ccthet_nucl(:,:,:,:,:)=0.0d0
1297 ddthet_nucl(:,:,:,:,:)=0.0d0
1298 eethet_nucl(:,:,:,:,:)=0.0d0
1299 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1300 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1305 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1306 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1307 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1308 read (ithep_nucl,*,end=111,err=111) &
1309 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1310 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1311 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1312 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1313 read (ithep_nucl,*,end=111,err=111) &
1314 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1315 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1316 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1321 !-------------------------------------------
1322 allocate(nlob(ntyp1)) !(ntyp1)
1323 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1324 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1325 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1332 gaussc(:,:,:,:)=0.0D0
1336 ! Read the parameters of the probability distribution/energy expression
1337 ! of the side chains.
1340 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1344 dsc_inv(i)=1.0D0/dsc(i)
1355 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1356 ((blower(k,l,1),l=1,k),k=1,3)
1357 censc(1,1,-i)=censc(1,1,i)
1358 censc(2,1,-i)=censc(2,1,i)
1359 censc(3,1,-i)=-censc(3,1,i)
1361 read (irotam,*,end=112,err=112) bsc(j,i)
1362 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1363 ((blower(k,l,j),l=1,k),k=1,3)
1364 censc(1,j,-i)=censc(1,j,i)
1365 censc(2,j,-i)=censc(2,j,i)
1366 censc(3,j,-i)=-censc(3,j,i)
1367 ! BSC is amplitude of Gaussian
1374 akl=akl+blower(k,m,j)*blower(l,m,j)
1378 if (((k.eq.3).and.(l.ne.3)) &
1379 .or.((l.eq.3).and.(k.ne.3))) then
1380 gaussc(k,l,j,-i)=-akl
1381 gaussc(l,k,j,-i)=-akl
1383 gaussc(k,l,j,-i)=akl
1384 gaussc(l,k,j,-i)=akl
1393 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1396 if (nlobi.gt.0) then
1398 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1399 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1400 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1401 'log h',(bsc(j,i),j=1,nlobi)
1402 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1403 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1405 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1406 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1409 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1410 write (iout,'(a,f10.4,4(16x,f10.4))') &
1411 'Center ',(bsc(j,i),j=1,nlobi)
1412 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1421 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1422 ! added by Urszula Kozlowska 07/11/2007
1424 !el Maximum number of SC local term fitting function coefficiants
1425 !el integer,parameter :: maxsccoef=65
1427 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1430 read (irotam,*,end=112,err=112)
1432 read (irotam,*,end=112,err=112)
1435 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1439 !---------reading nucleic acid parameters for rotamers-------------------
1440 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1441 do i=1,ntyp_molec(2)
1442 read (irotam_nucl,*,end=112,err=112)
1444 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1450 write (iout,*) "Base rotamer parameters"
1451 do i=1,ntyp_molec(2)
1452 write (iout,'(a)') restyp(i,2)
1453 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1458 ! Read the parameters of the probability distribution/energy expression
1459 ! of the side chains.
1461 write (2,*) "Start reading ROTAM_PDB"
1463 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1467 dsc_inv(i)=1.0D0/dsc(i)
1478 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1479 ((blower(k,l,1),l=1,k),k=1,3)
1481 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1482 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1483 ((blower(k,l,j),l=1,k),k=1,3)
1490 akl=akl+blower(k,m,j)*blower(l,m,j)
1500 write (2,*) "End reading ROTAM_PDB"
1506 ! Read torsional parameters in old format
1508 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1510 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1511 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1512 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1514 !el from energy module--------
1515 allocate(v1(nterm_old,ntortyp,ntortyp))
1516 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1517 !el---------------------------
1522 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1528 write (iout,'(/a/)') 'Torsional constants:'
1531 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1532 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1538 ! Read torsional parameters
1540 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1541 read (itorp,*,end=113,err=113) ntortyp
1542 !el from energy module---------
1543 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1544 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1546 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1547 allocate(vlor2(maxlor,ntortyp,ntortyp))
1548 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1549 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1551 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1552 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1553 !el---------------------------
1556 !el---------------------------
1558 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1560 itortyp(i)=-itortyp(-i)
1562 itortyp(ntyp1)=ntortyp
1563 itortyp(-ntyp1)=-ntortyp
1565 write (iout,*) 'ntortyp',ntortyp
1567 do j=-ntortyp+1,ntortyp-1
1568 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1570 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1571 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1574 do k=1,nterm(i,j,iblock)
1575 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1577 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1578 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1579 v0ij=v0ij+si*v1(k,i,j,iblock)
1581 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1582 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1583 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1585 do k=1,nlor(i,j,iblock)
1586 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1587 vlor2(k,i,j),vlor3(k,i,j)
1588 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1591 v0(-i,-j,iblock)=v0ij
1597 write (iout,'(/a/)') 'Torsional constants:'
1599 do i=-ntortyp,ntortyp
1600 do j=-ntortyp,ntortyp
1601 write (iout,*) 'ityp',i,' jtyp',j
1602 write (iout,*) 'Fourier constants'
1603 do k=1,nterm(i,j,iblock)
1604 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1607 write (iout,*) 'Lorenz constants'
1608 do k=1,nlor(i,j,iblock)
1609 write (iout,'(3(1pe15.5))') &
1610 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1616 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1618 ! 6/23/01 Read parameters for double torsionals
1620 !el from energy module------------
1621 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1622 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1623 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1624 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1625 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1626 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1627 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1628 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1629 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1630 !---------------------------------
1634 do j=-ntortyp+1,ntortyp-1
1635 do k=-ntortyp+1,ntortyp-1
1636 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1637 ! write (iout,*) "OK onelett",
1640 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1641 .or. t3.ne.toronelet(k)) then
1642 write (iout,*) "Error in double torsional parameter file",&
1645 call MPI_Finalize(Ierror)
1647 stop "Error in double torsional parameter file"
1649 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1650 ntermd_2(i,j,k,iblock)
1651 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1652 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1653 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1654 ntermd_1(i,j,k,iblock))
1655 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1656 ntermd_1(i,j,k,iblock))
1657 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1658 ntermd_1(i,j,k,iblock))
1659 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1660 ntermd_1(i,j,k,iblock))
1661 ! Martix of D parameters for one dimesional foureir series
1662 do l=1,ntermd_1(i,j,k,iblock)
1663 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1664 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1665 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1666 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1667 ! write(iout,*) "whcodze" ,
1668 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1670 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1671 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1672 v2s(m,l,i,j,k,iblock),&
1673 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1674 ! Martix of D parameters for two dimesional fourier series
1675 do l=1,ntermd_2(i,j,k,iblock)
1677 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1678 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1679 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1680 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1689 write (iout,*) 'Constants for double torsionals'
1692 do j=-ntortyp+1,ntortyp-1
1693 do k=-ntortyp+1,ntortyp-1
1694 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1695 ' nsingle',ntermd_1(i,j,k,iblock),&
1696 ' ndouble',ntermd_2(i,j,k,iblock)
1698 write (iout,*) 'Single angles:'
1699 do l=1,ntermd_1(i,j,k,iblock)
1700 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1701 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1702 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1703 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1706 write (iout,*) 'Pairs of angles:'
1707 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1708 do l=1,ntermd_2(i,j,k,iblock)
1709 write (iout,'(i5,20f10.5)') &
1710 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1713 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1714 do l=1,ntermd_2(i,j,k,iblock)
1715 write (iout,'(i5,20f10.5)') &
1716 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1717 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1726 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1727 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
1728 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
1729 !el from energy module---------
1730 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1731 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1733 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1734 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
1735 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
1736 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1738 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
1739 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1740 !el---------------------------
1743 !el--------------------
1744 read (itorp_nucl,*,end=113,err=113) &
1745 (itortyp_nucl(i),i=1,ntyp_molec(2))
1746 ! print *,itortyp_nucl(:)
1747 !c write (iout,*) 'ntortyp',ntortyp
1750 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
1751 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
1754 do k=1,nterm_nucl(i,j)
1755 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
1756 v0ij=v0ij+si*v1_nucl(k,i,j)
1759 do k=1,nlor_nucl(i,j)
1760 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
1761 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
1762 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
1768 ! Read of Side-chain backbone correlation parameters
1769 ! Modified 11 May 2012 by Adasko
1772 read (isccor,*,end=119,err=119) nsccortyp
1774 !el from module energy-------------
1775 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1776 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1777 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1778 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1779 !-----------------------------------
1781 !el from module energy-------------
1782 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1784 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1786 isccortyp(i)=-isccortyp(-i)
1788 iscprol=isccortyp(20)
1789 ! write (iout,*) 'ntortyp',ntortyp
1791 !c maxinter is maximum interaction sites
1792 !el from module energy---------
1793 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1794 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1795 -nsccortyp:nsccortyp))
1796 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1797 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1798 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1799 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1800 !-----------------------------------
1804 read (isccor,*,end=119,err=119) &
1805 nterm_sccor(i,j),nlor_sccor(i,j)
1811 nterm_sccor(-i,j)=nterm_sccor(i,j)
1812 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1813 nterm_sccor(i,-j)=nterm_sccor(i,j)
1814 do k=1,nterm_sccor(i,j)
1815 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1817 if (j.eq.iscprol) then
1818 if (i.eq.isccortyp(10)) then
1819 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1820 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1822 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1823 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1824 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1825 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1826 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1827 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1828 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1829 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1832 if (i.eq.isccortyp(10)) then
1833 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1834 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1836 if (j.eq.isccortyp(10)) then
1837 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1838 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1840 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1841 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1842 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1843 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1844 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1845 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1849 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1850 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1851 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1852 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1855 do k=1,nlor_sccor(i,j)
1856 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1857 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1858 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1859 (1+vlor3sccor(k,i,j)**2)
1861 v0sccor(l,i,j)=v0ijsccor
1862 v0sccor(l,-i,j)=v0ijsccor1
1863 v0sccor(l,i,-j)=v0ijsccor2
1864 v0sccor(l,-i,-j)=v0ijsccor3
1870 !el from module energy-------------
1871 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1873 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1874 ! write (iout,*) 'ntortyp',ntortyp
1876 !c maxinter is maximum interaction sites
1877 !el from module energy---------
1878 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1879 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1880 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1881 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1882 !-----------------------------------
1886 read (isccor,*,end=119,err=119) &
1887 nterm_sccor(i,j),nlor_sccor(i,j)
1891 do k=1,nterm_sccor(i,j)
1892 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1894 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1897 do k=1,nlor_sccor(i,j)
1898 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1899 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1900 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1901 (1+vlor3sccor(k,i,j)**2)
1903 v0sccor(l,i,j)=v0ijsccor !el ,iblock
1912 write (iout,'(/a/)') 'Torsional constants:'
1915 write (iout,*) 'ityp',i,' jtyp',j
1916 write (iout,*) 'Fourier constants'
1917 do k=1,nterm_sccor(i,j)
1918 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1920 write (iout,*) 'Lorenz constants'
1921 do k=1,nlor_sccor(i,j)
1922 write (iout,'(3(1pe15.5))') &
1923 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1930 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1931 ! interaction energy of the Gly, Ala, and Pro prototypes.
1936 write (iout,*) "Coefficients of the cumulants"
1938 read (ifourier,*) nloctyp
1939 !write(iout,*) "nloctyp",nloctyp
1940 !el from module energy-------
1941 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1942 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1943 allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1944 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1945 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1946 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1947 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1948 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1952 !--------------------------------
1955 read (ifourier,*,end=115,err=115)
1956 read (ifourier,*,end=115,err=115) (buse(ii),ii=1,13)
1958 write (iout,*) 'Type',i
1959 write (iout,'(a,i2,a,f10.5)') ('buse(',ii,')=',buse(ii),ii=1,13)
1967 B1tilde(1,i) = buse(3)
1968 B1tilde(2,i) =-buse(5)
1969 B1tilde(1,-i) =-buse(3)
1970 B1tilde(2,-i) =buse(5)
1971 ! buse1tilde(1,i)=0.0d0
1972 ! buse1tilde(2,i)=0.0d0
1992 Ctilde(1,1,i)=buse(7)
1993 Ctilde(1,2,i)=buse(9)
1994 Ctilde(2,1,i)=-buse(9)
1995 Ctilde(2,2,i)=buse(7)
1996 Ctilde(1,1,-i)=buse(7)
1997 Ctilde(1,2,-i)=-buse(9)
1998 Ctilde(2,1,-i)=buse(9)
1999 Ctilde(2,2,-i)=buse(7)
2001 ! Ctilde(1,1,i)=0.0d0
2002 ! Ctilde(1,2,i)=0.0d0
2003 ! Ctilde(2,1,i)=0.0d0
2004 ! Ctilde(2,2,i)=0.0d0
2017 Dtilde(1,1,i)=buse(6)
2018 Dtilde(1,2,i)=buse(8)
2019 Dtilde(2,1,i)=-buse(8)
2020 Dtilde(2,2,i)=buse(6)
2021 Dtilde(1,1,-i)=buse(6)
2022 Dtilde(1,2,-i)=-buse(8)
2023 Dtilde(2,1,-i)=buse(8)
2024 Dtilde(2,2,-i)=buse(6)
2026 ! Dtilde(1,1,i)=0.0d0
2027 ! Dtilde(1,2,i)=0.0d0
2028 ! Dtilde(2,1,i)=0.0d0
2029 ! Dtilde(2,2,i)=0.0d0
2030 EE(1,1,i)= buse(10)+buse(11)
2031 EE(2,2,i)=-buse(10)+buse(11)
2032 EE(2,1,i)= buse(12)-buse(13)
2033 EE(1,2,i)= buse(12)+buse(13)
2034 EE(1,1,-i)= buse(10)+buse(11)
2035 EE(2,2,-i)=-buse(10)+buse(11)
2036 EE(2,1,-i)=-buse(12)+buse(13)
2037 EE(1,2,-i)=-buse(12)-buse(13)
2043 ! ee(2,1,i)=ee(1,2,i)
2047 write (iout,*) 'Type',i
2049 write(iout,*) B1(1,i),B1(2,i)
2051 write(iout,*) B2(1,i),B2(2,i)
2054 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
2058 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
2062 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
2067 ! Read electrostatic-interaction parameters
2072 write (iout,'(/a)') 'Electrostatic interaction constants:'
2073 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2074 'IT','JT','APP','BPP','AEL6','AEL3'
2076 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2077 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2078 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2079 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2084 app (i,j)=epp(i,j)*rri*rri
2085 bpp (i,j)=-2.0D0*epp(i,j)*rri
2086 ael6(i,j)=elpp6(i,j)*4.2D0**6
2087 ael3(i,j)=elpp3(i,j)*4.2D0**3
2089 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2095 ! Read side-chain interaction parameters.
2097 !el from module energy - COMMON.INTERACT-------
2098 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2099 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2100 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2102 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2103 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2104 allocate(epslip(ntyp,ntyp))
2112 !--------------------------------
2114 read (isidep,*,end=117,err=117) ipot,expon
2115 if (ipot.lt.1 .or. ipot.gt.5) then
2116 write (iout,'(2a)') 'Error while reading SC interaction',&
2117 'potential file - unknown potential type.'
2119 call MPI_Finalize(Ierror)
2125 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2126 ', exponents are ',expon,2*expon
2127 ! goto (10,20,30,30,40) ipot
2129 !----------------------- LJ potential ---------------------------------
2131 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2132 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2133 (sigma0(i),i=1,ntyp)
2135 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2136 write (iout,'(a/)') 'The epsilon array:'
2137 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2138 write (iout,'(/a)') 'One-body parameters:'
2139 write (iout,'(a,4x,a)') 'residue','sigma'
2140 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2143 !----------------------- LJK potential --------------------------------
2145 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2146 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2147 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2149 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2150 write (iout,'(a/)') 'The epsilon array:'
2151 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2152 write (iout,'(/a)') 'One-body parameters:'
2153 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2154 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2158 !---------------------- GB or BP potential -----------------------------
2161 ! print *,"I AM in SCELE",scelemode
2162 if (scelemode.eq.0) then
2164 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2166 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2167 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2168 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2169 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2171 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2174 ! For the GB potential convert sigma'**2 into chi'
2177 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2181 write (iout,'(/a/)') 'Parameters of the BP potential:'
2182 write (iout,'(a/)') 'The epsilon array:'
2183 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2184 write (iout,'(/a)') 'One-body parameters:'
2185 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2187 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2188 chip(i),alp(i),i=1,ntyp)
2191 ! print *,ntyp,"NTYP"
2192 allocate(icharge(ntyp1))
2193 ! print *,ntyp,icharge(i)
2195 read (isidep,*) (icharge(i),i=1,ntyp)
2196 print *,ntyp,icharge(i)
2197 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2198 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2199 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2200 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2201 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2202 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2203 allocate(epsintab(ntyp,ntyp))
2204 allocate(dtail(2,ntyp,ntyp))
2205 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2206 allocate(wqdip(2,ntyp,ntyp))
2207 allocate(wstate(4,ntyp,ntyp))
2208 allocate(dhead(2,2,ntyp,ntyp))
2209 allocate(nstate(ntyp,ntyp))
2210 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2211 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2214 ! write (*,*) "Im in ALAB", i, " ", j
2216 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
2217 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
2218 chis(i,j),chis(j,i), &
2219 nstate(i,j),(wstate(k,i,j),k=1,4), &
2220 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
2221 dtail(1,i,j),dtail(2,i,j), &
2222 epshead(i,j),sig0head(i,j), &
2223 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
2224 alphapol(i,j),alphapol(j,i), &
2225 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
2226 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
2232 sigma(i,j) = sigma(j,i)
2233 sigmap1(i,j)=sigmap1(j,i)
2234 sigmap2(i,j)=sigmap2(j,i)
2235 sigiso1(i,j)=sigiso1(j,i)
2236 sigiso2(i,j)=sigiso2(j,i)
2237 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2238 nstate(i,j) = nstate(j,i)
2239 dtail(1,i,j) = dtail(1,j,i)
2240 dtail(2,i,j) = dtail(2,j,i)
2242 alphasur(k,i,j) = alphasur(k,j,i)
2243 wstate(k,i,j) = wstate(k,j,i)
2244 alphiso(k,i,j) = alphiso(k,j,i)
2247 dhead(2,1,i,j) = dhead(1,1,j,i)
2248 dhead(2,2,i,j) = dhead(1,2,j,i)
2249 dhead(1,1,i,j) = dhead(2,1,j,i)
2250 dhead(1,2,i,j) = dhead(2,2,j,i)
2252 epshead(i,j) = epshead(j,i)
2253 sig0head(i,j) = sig0head(j,i)
2256 wqdip(k,i,j) = wqdip(k,j,i)
2259 wquad(i,j) = wquad(j,i)
2260 epsintab(i,j) = epsintab(j,i)
2261 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2266 !--------------------- GBV potential -----------------------------------
2268 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2269 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2270 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2271 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2273 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2274 write (iout,'(a/)') 'The epsilon array:'
2275 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2276 write (iout,'(/a)') 'One-body parameters:'
2277 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2278 's||/s_|_^2',' chip ',' alph '
2279 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2280 sigii(i),chip(i),alp(i),i=1,ntyp)
2283 write(iout,*)"Wrong ipot"
2289 !-----------------------------------------------------------------------
2290 ! Calculate the "working" parameters of SC interactions.
2292 !el from module energy - COMMON.INTERACT-------
2293 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2294 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2295 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2296 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2297 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2298 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2300 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2306 if (scelemode.eq.0) then
2315 sc_aa_tube_par(:)=0.0d0
2316 sc_bb_tube_par(:)=0.0d0
2318 !--------------------------------
2323 epslip(i,j)=epslip(j,i)
2326 if (scelemode.eq.0) then
2329 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2330 sigma(j,i)=sigma(i,j)
2331 rs0(i,j)=dwa16*sigma(i,j)
2336 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2337 'Working parameters of the SC interactions:',&
2338 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2343 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2345 ! print *,"SIGMA ZLA?",sigma(i,j)
2353 sigeps=dsign(1.0D0,epsij)
2355 aa_aq(i,j)=epsij*rrij*rrij
2356 print *,"ADASKO",epsij,rrij,expon
2357 bb_aq(i,j)=-sigeps*epsij*rrij
2358 aa_aq(j,i)=aa_aq(i,j)
2359 bb_aq(j,i)=bb_aq(i,j)
2360 epsijlip=epslip(i,j)
2361 sigeps=dsign(1.0D0,epsijlip)
2362 epsijlip=dabs(epsijlip)
2363 aa_lip(i,j)=epsijlip*rrij*rrij
2364 bb_lip(i,j)=-sigeps*epsijlip*rrij
2365 aa_lip(j,i)=aa_lip(i,j)
2366 bb_lip(j,i)=bb_lip(i,j)
2367 !C write(iout,*) aa_lip
2368 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2369 sigt1sq=sigma0(i)**2
2370 sigt2sq=sigma0(j)**2
2373 ratsig1=sigt2sq/sigt1sq
2374 ratsig2=1.0D0/ratsig1
2375 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2376 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2377 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2381 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2382 sigmaii(i,j)=rsum_max
2383 sigmaii(j,i)=rsum_max
2385 ! sigmaii(i,j)=r0(i,j)
2386 ! sigmaii(j,i)=r0(i,j)
2388 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2389 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2390 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2391 augm(i,j)=epsij*r_augm**(2*expon)
2392 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2399 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2400 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2401 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2406 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2407 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2408 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2409 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2410 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2411 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2412 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2413 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2414 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2415 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2416 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2417 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2418 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2419 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2420 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2421 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2430 read (isidep_nucl,*) ipot_nucl
2431 ! print *,"TU?!",ipot_nucl
2432 if (ipot_nucl.eq.1) then
2433 do i=1,ntyp_molec(2)
2434 do j=i,ntyp_molec(2)
2435 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2436 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2440 do i=1,ntyp_molec(2)
2441 do j=i,ntyp_molec(2)
2442 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2443 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2444 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2448 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2449 do i=1,ntyp_molec(2)
2450 do j=i,ntyp_molec(2)
2451 rrij=sigma_nucl(i,j)
2455 epsij=4*eps_nucl(i,j)
2456 sigeps=dsign(1.0D0,epsij)
2458 aa_nucl(i,j)=epsij*rrij*rrij
2459 bb_nucl(i,j)=-sigeps*epsij*rrij
2460 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2461 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2462 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2463 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2464 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2465 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2468 aa_nucl(i,j)=aa_nucl(j,i)
2469 bb_nucl(i,j)=bb_nucl(j,i)
2470 ael3_nucl(i,j)=ael3_nucl(j,i)
2471 ael6_nucl(i,j)=ael6_nucl(j,i)
2472 ael63_nucl(i,j)=ael63_nucl(j,i)
2473 ael32_nucl(i,j)=ael32_nucl(j,i)
2474 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2475 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2476 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2477 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2478 eps_nucl(i,j)=eps_nucl(j,i)
2479 sigma_nucl(i,j)=sigma_nucl(j,i)
2480 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2484 write(iout,*) "tube param"
2485 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2486 ccavtubpep,dcavtubpep,tubetranenepep
2487 sigmapeptube=sigmapeptube**6
2488 sigeps=dsign(1.0D0,epspeptube)
2489 epspeptube=dabs(epspeptube)
2490 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2491 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2492 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2494 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2495 ccavtub(i),dcavtub(i),tubetranene(i)
2496 sigmasctube=sigmasctube**6
2497 sigeps=dsign(1.0D0,epssctube)
2498 epssctube=dabs(epssctube)
2499 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2500 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2501 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2503 !-----------------READING SC BASE POTENTIALS-----------------------------
2504 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2505 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2506 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2507 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2508 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2509 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2510 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2511 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2512 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2513 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2514 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2515 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2516 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2517 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2518 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2519 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2522 do i=1,ntyp_molec(1)
2523 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2524 write (*,*) "Im in ", i, " ", j
2525 read(isidep_scbase,*) &
2526 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2527 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2528 write(*,*) "eps",eps_scbase(i,j)
2529 read(isidep_scbase,*) &
2530 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2531 chis_scbase(i,j,1),chis_scbase(i,j,2)
2532 read(isidep_scbase,*) &
2533 dhead_scbasei(i,j), &
2534 dhead_scbasej(i,j), &
2535 rborn_scbasei(i,j),rborn_scbasej(i,j)
2536 read(isidep_scbase,*) &
2537 (wdipdip_scbase(k,i,j),k=1,3), &
2538 (wqdip_scbase(k,i,j),k=1,2)
2539 read(isidep_scbase,*) &
2540 alphapol_scbase(i,j), &
2541 epsintab_scbase(i,j)
2544 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2545 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2547 do i=1,ntyp_molec(1)
2548 do j=1,ntyp_molec(2)-1
2549 epsij=eps_scbase(i,j)
2550 rrij=sigma_scbase(i,j)
2555 sigeps=dsign(1.0D0,epsij)
2557 aa_scbase(i,j)=epsij*rrij*rrij
2558 bb_scbase(i,j)=-sigeps*epsij*rrij
2561 !-----------------READING PEP BASE POTENTIALS-------------------
2562 allocate(eps_pepbase(ntyp_molec(2)))
2563 allocate(sigma_pepbase(ntyp_molec(2)))
2564 allocate(chi_pepbase(ntyp_molec(2),2))
2565 allocate(chipp_pepbase(ntyp_molec(2),2))
2566 allocate(alphasur_pepbase(4,ntyp_molec(2)))
2567 allocate(sigmap1_pepbase(ntyp_molec(2)))
2568 allocate(sigmap2_pepbase(ntyp_molec(2)))
2569 allocate(chis_pepbase(ntyp_molec(2),2))
2570 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
2573 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2574 write (*,*) "Im in ", i, " ", j
2575 read(isidep_pepbase,*) &
2576 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
2577 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
2578 write(*,*) "eps",eps_pepbase(j)
2579 read(isidep_pepbase,*) &
2580 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
2581 chis_pepbase(j,1),chis_pepbase(j,2)
2582 read(isidep_pepbase,*) &
2583 (wdipdip_pepbase(k,j),k=1,3)
2585 allocate(aa_pepbase(ntyp_molec(2)))
2586 allocate(bb_pepbase(ntyp_molec(2)))
2588 do j=1,ntyp_molec(2)-1
2589 epsij=eps_pepbase(j)
2590 rrij=sigma_pepbase(j)
2595 sigeps=dsign(1.0D0,epsij)
2597 aa_pepbase(j)=epsij*rrij*rrij
2598 bb_pepbase(j)=-sigeps*epsij*rrij
2600 !--------------READING SC PHOSPHATE-------------------------------------
2601 allocate(eps_scpho(ntyp_molec(1)))
2602 allocate(sigma_scpho(ntyp_molec(1)))
2603 allocate(chi_scpho(ntyp_molec(1),2))
2604 allocate(chipp_scpho(ntyp_molec(1),2))
2605 allocate(alphasur_scpho(4,ntyp_molec(1)))
2606 allocate(sigmap1_scpho(ntyp_molec(1)))
2607 allocate(sigmap2_scpho(ntyp_molec(1)))
2608 allocate(chis_scpho(ntyp_molec(1),2))
2609 allocate(wqq_scpho(ntyp_molec(1)))
2610 allocate(wqdip_scpho(2,ntyp_molec(1)))
2611 allocate(alphapol_scpho(ntyp_molec(1)))
2612 allocate(epsintab_scpho(ntyp_molec(1)))
2613 allocate(dhead_scphoi(ntyp_molec(1)))
2614 allocate(rborn_scphoi(ntyp_molec(1)))
2615 allocate(rborn_scphoj(ntyp_molec(1)))
2616 allocate(alphi_scpho(ntyp_molec(1)))
2620 do j=1,ntyp_molec(1) ! without U then we will take T for U
2621 write (*,*) "Im in scpho ", i, " ", j
2622 read(isidep_scpho,*) &
2623 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
2624 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
2625 write(*,*) "eps",eps_scpho(j)
2626 read(isidep_scpho,*) &
2627 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
2628 chis_scpho(j,1),chis_scpho(j,2)
2629 read(isidep_scpho,*) &
2630 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
2631 read(isidep_scpho,*) &
2632 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
2636 allocate(aa_scpho(ntyp_molec(1)))
2637 allocate(bb_scpho(ntyp_molec(1)))
2639 do j=1,ntyp_molec(1)
2646 sigeps=dsign(1.0D0,epsij)
2648 aa_scpho(j)=epsij*rrij*rrij
2649 bb_scpho(j)=-sigeps*epsij*rrij
2653 read(isidep_peppho,*) &
2654 eps_peppho,sigma_peppho
2655 read(isidep_peppho,*) &
2656 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
2657 read(isidep_peppho,*) &
2658 (wqdip_peppho(k),k=1,2)
2666 sigeps=dsign(1.0D0,epsij)
2668 aa_peppho=epsij*rrij*rrij
2669 bb_peppho=-sigeps*epsij*rrij
2672 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2677 ! Define the SC-p interaction constants (hard-coded; old style)
2680 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2682 ! aad(i,1)=0.3D0*4.0D0**12
2683 ! Following line for constants currently implemented
2684 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2685 aad(i,1)=1.5D0*4.0D0**12
2686 ! aad(i,1)=0.17D0*5.6D0**12
2688 ! "Soft" SC-p repulsion
2690 ! Following line for constants currently implemented
2691 ! aad(i,1)=0.3D0*4.0D0**6
2692 ! "Hard" SC-p repulsion
2693 bad(i,1)=3.0D0*4.0D0**6
2694 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2703 ! 8/9/01 Read the SC-p interaction constants from file
2706 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2709 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2710 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2711 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2712 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2716 write (iout,*) "Parameters of SC-p interactions:"
2718 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2719 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2724 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2726 do i=1,ntyp_molec(2)
2727 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
2729 do i=1,ntyp_molec(2)
2730 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2731 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2733 r0pp=1.12246204830937298142*5.16158
2739 ! Define the constants of the disulfide bridge
2743 ! Old arbitrary potential - commented out.
2748 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2749 ! energy surface of diethyl disulfide.
2750 ! A. Liwo and U. Kozlowska, 11/24/03
2767 write (iout,'(/a)') "Disulfide bridge parameters:"
2768 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2769 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2770 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2771 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2774 if (shield_mode.gt.0) then
2775 pi=4.0D0*datan(1.0D0)
2776 !C VSolvSphere the volume of solving sphere
2778 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
2779 !C there will be no distinction between proline peptide group and normal peptide
2780 !C group in case of shielding parameters
2781 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
2782 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
2783 write (iout,*) VSolvSphere,VSolvSphere_div
2784 !C long axis of side chain
2786 long_r_sidechain(i)=vbldsc0(1,i)
2787 ! if (scelemode.eq.0) then
2788 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
2789 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
2791 ! short_r_sidechain(i)=sigma(i,i)
2793 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
2800 111 write (iout,*) "Error reading bending energy parameters."
2802 112 write (iout,*) "Error reading rotamer energy parameters."
2804 113 write (iout,*) "Error reading torsional energy parameters."
2806 114 write (iout,*) "Error reading double torsional energy parameters."
2808 115 write (iout,*) &
2809 "Error reading cumulant (multibody energy) parameters."
2811 116 write (iout,*) "Error reading electrostatic energy parameters."
2813 117 write (iout,*) "Error reading side chain interaction parameters."
2815 118 write (iout,*) "Error reading SCp interaction parameters."
2817 119 write (iout,*) "Error reading SCCOR parameters"
2820 call MPI_Finalize(Ierror)
2824 end subroutine parmread
2826 !-----------------------------------------------------------------------------
2828 !-----------------------------------------------------------------------------
2829 subroutine printmat(ldim,m,n,iout,key,a)
2832 character(len=3),dimension(n) :: key
2833 real(kind=8),dimension(ldim,n) :: a
2835 integer :: i,j,k,m,iout,nlim
2839 write (iout,1000) (key(k),k=i,nlim)
2841 1000 format (/5x,8(6x,a3))
2842 1020 format (/80(1h-)/)
2844 write (iout,1010) key(j),(a(j,k),k=i,nlim)
2847 1010 format (a3,2x,8(f9.4))
2849 end subroutine printmat
2850 !-----------------------------------------------------------------------------
2852 !-----------------------------------------------------------------------------
2854 ! Read the PDB file and convert the peptide geometry into virtual-chain
2857 use energy_data, only: itype,istype
2861 use control, only: rescode,sugarcode
2862 ! implicit real*8 (a-h,o-z)
2863 ! include 'DIMENSIONS'
2864 ! include 'COMMON.LOCAL'
2865 ! include 'COMMON.VAR'
2866 ! include 'COMMON.CHAIN'
2867 ! include 'COMMON.INTERACT'
2868 ! include 'COMMON.IOUNITS'
2869 ! include 'COMMON.GEO'
2870 ! include 'COMMON.NAMES'
2871 ! include 'COMMON.CONTROL'
2872 ! include 'COMMON.DISTFIT'
2873 ! include 'COMMON.SETUP'
2874 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
2876 logical :: lprn=.true.,fail
2877 real(kind=8),dimension(3) :: e1,e2,e3
2878 real(kind=8) :: dcj,efree_temp
2879 character(len=3) :: seq,res
2880 character(len=5) :: atom
2881 character(len=80) :: card
2882 real(kind=8),dimension(3,20) :: sccor
2883 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
2884 integer :: isugar,molecprev,firstion
2885 character*1 :: sugar
2887 real(kind=8),dimension(3) :: ccc
2889 integer,dimension(2,maxres/3) :: hfrag_alloc
2890 integer,dimension(4,maxres/3) :: bfrag_alloc
2891 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2892 real(kind=8),dimension(:,:), allocatable :: c_temporary
2893 integer,dimension(:,:) , allocatable :: itype_temporary
2894 integer,dimension(:),allocatable :: istype_temp
2901 ! write (2,*) "UNRES_PDB",unres_pdb
2921 !-----------------------------
2922 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2923 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2924 if(.not. allocated(istype)) allocate(istype(maxres))
2926 read (ipdbin,'(a80)',end=10) card
2927 ! write (iout,'(a)') card
2928 if (card(:5).eq.'HELIX') then
2931 read(card(22:25),*) hfrag(1,nhfrag)
2932 read(card(34:37),*) hfrag(2,nhfrag)
2934 if (card(:5).eq.'SHEET') then
2937 read(card(24:26),*) bfrag(1,nbfrag)
2938 read(card(35:37),*) bfrag(2,nbfrag)
2939 !rc----------------------------------------
2940 !rc to be corrected !!!
2941 bfrag(3,nbfrag)=bfrag(1,nbfrag)
2942 bfrag(4,nbfrag)=bfrag(2,nbfrag)
2943 !rc----------------------------------------
2945 if (card(:3).eq.'END') then
2947 else if (card(:3).eq.'TER') then
2952 itype(ires_old,molecule)=ntyp1_molec(molecule)
2953 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
2954 nres_molec(molecule)=nres_molec(molecule)+2
2956 ! write (iout,*) "Chain ended",ires,ishift,ires_old
2959 dc(j,ires)=sccor(j,iii)
2962 call sccenter(ires,iii,sccor)
2968 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2969 ! Fish out the ATOM cards.
2970 ! write(iout,*) 'card',card(1:20)
2971 if (index(card(1:4),'ATOM').gt.0) then
2972 read (card(12:16),*) atom
2973 ! write (iout,*) "! ",atom," !",ires
2974 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
2975 read (card(23:26),*) ires
2976 read (card(18:20),'(a3)') res
2977 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
2978 ! & " ires_old",ires_old
2979 ! write (iout,*) "ishift",ishift," ishift1",ishift1
2980 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2981 if (ires-ishift+ishift1.ne.ires_old) then
2982 ! Calculate the CM of the preceding residue.
2983 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2985 ! write (iout,*) "Calculating sidechain center iii",iii
2988 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2991 call sccenter(ires_old,iii,sccor)
2995 ! Start new residue.
2996 if (res.eq.'Cl-' .or. res.eq.'Na+') then
2999 else if (ibeg.eq.1) then
3000 write (iout,*) "BEG ires",ires
3002 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3005 nres_molec(molecule)=nres_molec(molecule)+1
3007 ires=ires-ishift+ishift1
3009 ! write (iout,*) "ishift",ishift," ires",ires,&
3010 ! " ires_old",ires_old
3012 else if (ibeg.eq.2) then
3014 ishift=-ires_old+ires-1 !!!!!
3015 ishift1=ishift1-1 !!!!!
3016 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3017 ires=ires-ishift+ishift1
3018 ! print *,ires,ishift,ishift1
3022 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
3023 ires=ires-ishift+ishift1
3026 ! print *,'atom',ires,atom
3027 if (res.eq.'ACE' .or. res.eq.'NHE') then
3030 if (atom.eq.'CA '.or.atom.eq.'N ') then
3032 itype(ires,molecule)=rescode(ires,res,0,molecule)
3034 ! nres_molec(molecule)=nres_molec(molecule)+1
3037 itype(ires,molecule)=rescode(ires,res(2:3),0,molecule)
3038 ! nres_molec(molecule)=nres_molec(molecule)+1
3039 read (card(19:19),'(a1)') sugar
3040 isugar=sugarcode(sugar,ires)
3041 ! if (ibeg.eq.1) then
3045 ! print *,"ires=",ires,istype(ires)
3051 ires=ires-ishift+ishift1
3053 ! write (iout,*) "ires_old",ires_old," ires",ires
3054 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3057 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3058 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
3059 res.eq.'NHE'.and.atom(:2).eq.'HN') then
3060 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3061 ! print *,ires,ishift,ishift1
3062 ! write (iout,*) "backbone ",atom
3064 write (iout,'(2i3,2x,a,3f8.3)') &
3065 ires,itype(ires,1),res,(c(j,ires),j=1,3)
3068 nres_molec(molecule)=nres_molec(molecule)+1
3070 sccor(j,iii)=c(j,ires)
3072 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
3073 atom.eq."C2'" .or. atom.eq."C3'" &
3074 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
3075 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
3076 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
3077 ! print *,ires,ishift,ishift1
3081 ! sccor(j,iii)=c(j,ires)
3084 c(j,ires)=c(j,ires)+ccc(j)/5.0
3086 print *,counter,molecule
3087 if (counter.eq.5) then
3089 nres_molec(molecule)=nres_molec(molecule)+1
3092 ! sccor(j,iii)=c(j,ires)
3096 ! print *, "ATOM",atom(1:3)
3097 else if (atom.eq."C5'") then
3098 read (card(19:19),'(a1)') sugar
3099 isugar=sugarcode(sugar,ires)
3104 ! print *,ires,istype(ires)
3108 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
3109 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3110 nres_molec(molecule)=nres_molec(molecule)+1
3111 print *,"nres_molec(molecule)",nres_molec(molecule),ires
3115 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3117 else if ((atom.eq."C1'").and.unres_pdb) then
3119 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3120 ! write (*,*) card(23:27),ires,itype(ires,1)
3121 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
3122 atom.ne.'N' .and. atom.ne.'C' .and. &
3123 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
3124 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
3125 .and. atom.ne.'P '.and. &
3126 atom(1:1).ne.'H' .and. &
3127 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
3129 ! write (iout,*) "sidechain ",atom
3130 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3131 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
3132 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3134 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3137 else if ((ions).and.(card(1:6).eq.'HETATM')) then
3138 if (firstion.eq.0) then
3142 dc(j,ires)=sccor(j,iii)
3145 call sccenter(ires,iii,sccor)
3148 read (card(12:16),*) atom
3149 print *,"HETATOM", atom
3150 read (card(18:20),'(a3)') res
3151 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3152 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3153 .or.(atom(1:2).eq.'K ')) &
3156 if (molecule.ne.5) molecprev=molecule
3158 nres_molec(molecule)=nres_molec(molecule)+1
3159 print *,"HERE",nres_molec(molecule)
3161 itype(ires,molecule)=rescode(ires,res,0,molecule)
3162 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3166 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3167 if (ires.eq.0) return
3168 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3171 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
3173 nres_molec(molecule)=nres_molec(molecule)-2
3174 print *,'I have',nres, nres_molec(:)
3176 do k=1,4 ! ions are without dummy
3177 if (nres_molec(k).eq.0) cycle
3179 ! write (iout,*) i,itype(i,1)
3180 ! if (itype(i,1).eq.ntyp1) then
3181 ! write (iout,*) "dummy",i,itype(i,1)
3183 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3184 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3188 if (itype(i,k).eq.ntyp1_molec(k)) then
3189 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3190 if (itype(i+2,k).eq.0) then
3191 ! print *,"masz sieczke"
3193 if (itype(i+2,j).ne.0) then
3195 itype(i+1,j)=ntyp1_molec(j)
3196 nres_molec(k)=nres_molec(k)-1
3197 nres_molec(j)=nres_molec(j)+1
3203 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3204 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3205 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3206 ! if (unres_pdb) then
3207 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3208 ! print *,i,'tu dochodze'
3209 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3217 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3221 dcj=(c(j,i-2)-c(j,i-3))/2.0
3222 if (dcj.eq.0) dcj=1.23591524223
3227 else !itype(i+1,1).eq.ntyp1
3228 ! if (unres_pdb) then
3229 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3230 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3237 c(j,i)=c(j,i+1)-1.9d0*e2(j)
3241 dcj=(c(j,i+3)-c(j,i+2))/2.0
3242 if (dcj.eq.0) dcj=1.23591524223
3247 endif !itype(i+1,1).eq.ntyp1
3248 endif !itype.eq.ntyp1
3252 ! Calculate the CM of the last side chain.
3256 dc(j,ires)=sccor(j,iii)
3259 call sccenter(ires,iii,sccor)
3265 ! print *,"molecule",molecule
3266 if ((itype(nres,1).ne.10)) then
3268 if (molecule.eq.5) molecule=molecprev
3269 itype(nres,molecule)=ntyp1_molec(molecule)
3270 nres_molec(molecule)=nres_molec(molecule)+1
3271 ! if (unres_pdb) then
3272 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3273 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3280 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3284 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3285 c(j,nres)=c(j,nres-1)+dcj
3286 c(j,2*nres)=c(j,nres)
3290 ! print *,'I have',nres, nres_molec(:)
3292 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3294 if (nres.ne.nres0) then
3295 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3297 stop "Error nres value in WHAM input"
3300 !---------------------------------
3301 !el reallocate tables
3304 ! hfrag_alloc(j,i)=hfrag(j,i)
3307 ! bfrag_alloc(j,i)=bfrag(j,i)
3313 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3314 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3315 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3316 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3320 ! hfrag(j,i)=hfrag_alloc(j,i)
3325 ! bfrag(j,i)=bfrag_alloc(j,i)
3328 !el end reallocate tables
3329 !---------------------------------
3337 c(j,2*nres)=c(j,nres)
3340 if (itype(1,1).eq.ntyp1) then
3344 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3345 call refsys(2,3,4,e1,e2,e3,fail)
3352 c(j,1)=c(j,2)-1.9d0*e2(j)
3356 dcj=(c(j,4)-c(j,3))/2.0
3362 ! First lets assign correct dummy to correct type of chain
3364 if (itype(1,1).eq.ntyp1) then
3365 if (itype(2,1).eq.0) then
3367 if (itype(2,j).ne.0) then
3369 itype(1,j)=ntyp1_molec(j)
3370 nres_molec(1)=nres_molec(1)-1
3371 nres_molec(j)=nres_molec(j)+1
3378 print *,'I have',nres, nres_molec(:)
3380 ! Copy the coordinates to reference coordinates
3386 ! Calculate internal coordinates.
3388 write (iout,'(/a)') &
3389 "Cartesian coordinates of the reference structure"
3390 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3391 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3393 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3394 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3395 (c(j,ires+nres),j=1,3)
3398 ! znamy już nres wiec mozna alokowac tablice
3399 ! Calculate internal coordinates.
3400 if(me.eq.king.or..not.out1file)then
3401 write (iout,'(a)') &
3402 "Backbone and SC coordinates as read from the PDB"
3404 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
3405 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3406 (c(j,nres+ires),j=1,3)
3409 ! NOW LETS ROCK! SORTING
3410 allocate(c_temporary(3,2*nres))
3411 allocate(itype_temporary(nres,5))
3412 if (.not.allocated(molnum)) allocate(molnum(nres+1))
3413 if (.not.allocated(istype)) write(iout,*) &
3414 "SOMETHING WRONG WITH ISTYTPE"
3415 allocate(istype_temp(nres))
3416 itype_temporary(:,:)=0
3420 if (itype(i,k).ne.0) then
3422 c_temporary(j,seqalingbegin)=c(j,i)
3423 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3426 itype_temporary(seqalingbegin,k)=itype(i,k)
3427 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3428 istype_temp(seqalingbegin)=istype(i)
3429 molnum(seqalingbegin)=k
3430 seqalingbegin=seqalingbegin+1
3436 c(j,i)=c_temporary(j,i)
3441 itype(i,k)=itype_temporary(i,k)
3442 istype(i)=istype_temp(i)
3445 ! if (itype(1,1).eq.ntyp1) then
3448 ! if (unres_pdb) then
3449 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3450 ! call refsys(2,3,4,e1,e2,e3,fail)
3457 ! c(j,1)=c(j,2)-1.9d0*e2(j)
3461 ! dcj=(c(j,4)-c(j,3))/2.0
3463 ! c(j,nres+1)=c(j,1)
3469 write (iout,'(/a)') &
3470 "Cartesian coordinates of the reference structure after sorting"
3471 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3472 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3474 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3475 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3476 (c(j,ires+nres),j=1,3)
3480 ! print *,seqalingbegin,nres
3481 if(.not.allocated(vbld)) then
3482 allocate(vbld(2*nres))
3487 if(.not.allocated(vbld_inv)) then
3488 allocate(vbld_inv(2*nres))
3494 if(.not.allocated(theta)) then
3495 allocate(theta(nres+2))
3499 if(.not.allocated(phi)) allocate(phi(nres+2))
3500 if(.not.allocated(alph)) allocate(alph(nres+2))
3501 if(.not.allocated(omeg)) allocate(omeg(nres+2))
3502 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3503 if(.not.allocated(phiref)) allocate(phiref(nres+2))
3504 if(.not.allocated(costtab)) allocate(costtab(nres))
3505 if(.not.allocated(sinttab)) allocate(sinttab(nres))
3506 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3507 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3508 if(.not.allocated(xxref)) allocate(xxref(nres))
3509 if(.not.allocated(yyref)) allocate(yyref(nres))
3510 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3511 if(.not.allocated(dc_norm)) then
3512 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3513 allocate(dc_norm(3,0:2*nres+2))
3517 call int_from_cart(.true.,.false.)
3518 call sc_loc_geom(.false.)
3520 thetaref(i)=theta(i)
3530 dc(j,i)=c(j,i+1)-c(j,i)
3531 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3536 dc(j,i+nres)=c(j,i+nres)-c(j,i)
3537 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3539 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3543 ! Copy the coordinates to reference coordinates
3544 ! Splits to single chain if occurs
3548 ! cref(j,i,cou)=c(j,i)
3552 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
3553 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
3554 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3555 !-----------------------------
3559 write (iout,*) "symetr", symetr
3562 !c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3564 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
3567 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3572 cref(j,i,cou)=c(j,i)
3573 cref(j,i+nres,cou)=c(j,i+nres)
3575 chain_rep(j,lll,kkk)=c(j,i)
3576 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
3580 write (iout,*) chain_length
3581 if (chain_length.eq.0) chain_length=nres
3583 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
3584 chain_rep(j,chain_length+nres,symetr) &
3585 =chain_rep(j,chain_length+nres,1)
3588 ! write (iout,*) "spraw lancuchy",chain_length,symetr
3590 ! do kkk=1,chain_length
3591 ! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
3595 ! makes copy of chains
3596 write (iout,*) "symetr", symetr
3601 if (symetr.gt.1) then
3608 write(iout,*) (tabperm(i,kkk),kkk=1,4)
3614 ! write (iout,*) i,icha
3615 do lll=1,chain_length
3617 if (cou.le.nres) then
3619 kupa=mod(lll,chain_length)
3620 iprzes=(kkk-1)*chain_length+lll
3621 if (kupa.eq.0) kupa=chain_length
3622 ! write (iout,*) "kupa", kupa
3623 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
3624 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
3631 !-koniec robienia kopii
3634 write (iout,*) "nowa struktura", nperm
3636 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
3638 cref(3,i,kkk),cref(1,nres+i,kkk),&
3639 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
3641 100 format (//' alpha-carbon coordinates ',&
3642 ' centroid coordinates'/ &
3643 ' ', 6X,'X',11X,'Y',11X,'Z', &
3644 10X,'X',11X,'Y',11X,'Z')
3645 110 format (a,'(',i3,')',6f12.5)
3651 bfrag(i,j)=bfrag(i,j)-ishift
3657 hfrag(i,j)=hfrag(i,j)-ishift
3663 end subroutine readpdb
3664 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3665 !-----------------------------------------------------------------------------
3667 !-----------------------------------------------------------------------------
3668 subroutine read_control
3682 use random, only: random_init
3683 ! implicit real*8 (a-h,o-z)
3684 ! include 'DIMENSIONS'
3686 use prng, only:prng_restart
3688 logical :: OKRandom!, prng_restart
3691 ! include 'COMMON.IOUNITS'
3692 ! include 'COMMON.TIME1'
3693 ! include 'COMMON.THREAD'
3694 ! include 'COMMON.SBRIDGE'
3695 ! include 'COMMON.CONTROL'
3696 ! include 'COMMON.MCM'
3697 ! include 'COMMON.MAP'
3698 ! include 'COMMON.HEADER'
3699 ! include 'COMMON.CSA'
3700 ! include 'COMMON.CHAIN'
3701 ! include 'COMMON.MUCA'
3702 ! include 'COMMON.MD'
3703 ! include 'COMMON.FFIELD'
3704 ! include 'COMMON.INTERACT'
3705 ! include 'COMMON.SETUP'
3706 !el integer :: KDIAG,ICORFL,IXDR
3707 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
3708 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
3709 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
3710 ! character(len=80) :: ucase
3711 character(len=640) :: controlcard
3713 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
3719 read (INP,'(a)') titel
3720 call card_concat(controlcard,.true.)
3721 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
3722 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
3723 call reada(controlcard,'SEED',seed,0.0D0)
3724 call random_init(seed)
3725 ! Set up the time limit (caution! The time must be input in minutes!)
3726 read_cart=index(controlcard,'READ_CART').gt.0
3727 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3728 call readi(controlcard,'SYM',symetr,1)
3729 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
3730 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
3731 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
3732 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
3733 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
3734 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
3735 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
3736 call reada(controlcard,'DRMS',drms,0.1D0)
3737 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3738 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
3739 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
3740 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
3741 write (iout,'(a,f10.1)')'DRMS = ',drms
3742 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
3743 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
3745 call readi(controlcard,'NZ_START',nz_start,0)
3746 call readi(controlcard,'NZ_END',nz_end,0)
3747 call readi(controlcard,'IZ_SC',iz_sc,0)
3748 timlim=60.0D0*timlim
3749 safety = 60.0d0*safety
3752 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3753 !C SHIELD keyword sets if the shielding effect of side-chains is used
3754 !C 0 denots no shielding is used all peptide are equally despite the
3755 !C solvent accesible area
3756 !C 1 the newly introduced function
3757 !C 2 reseved for further possible developement
3758 call readi(controlcard,'SHIELD',shield_mode,0)
3759 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3760 write(iout,*) "shield_mode",shield_mode
3761 !C Varibles set size of box
3762 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
3763 protein=index(controlcard,"PROTEIN").gt.0
3764 ions=index(controlcard,"IONS").gt.0
3765 nucleic=index(controlcard,"NUCLEIC").gt.0
3766 write (iout,*) "with_theta_constr ",with_theta_constr
3767 AFMlog=(index(controlcard,'AFM'))
3768 selfguide=(index(controlcard,'SELFGUIDE'))
3769 print *,'AFMlog',AFMlog,selfguide,"KUPA"
3770 call readi(controlcard,'GENCONSTR',genconstr,0)
3771 call reada(controlcard,'BOXX',boxxsize,100.0d0)
3772 call reada(controlcard,'BOXY',boxysize,100.0d0)
3773 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3774 call readi(controlcard,'TUBEMOD',tubemode,0)
3775 print *,"SCELE",scelemode
3776 call readi(controlcard,"SCELEMODE",scelemode,0)
3777 print *,"SCELE",scelemode
3779 ! elemode = 0 is orignal UNRES electrostatics
3780 ! elemode = 1 is "Momo" potentials in progress
3781 ! elemode = 2 is in development EVALD
3782 write (iout,*) TUBEmode,"TUBEMODE"
3783 if (TUBEmode.gt.0) then
3784 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3785 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3786 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3787 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3788 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3789 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3790 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3791 buftubebot=bordtubebot+tubebufthick
3792 buftubetop=bordtubetop-tubebufthick
3795 ! CUTOFFF ON ELECTROSTATICS
3796 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
3797 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3798 write(iout,*) "R_CUT_ELE=",r_cut_ele
3799 ! Lipidic parameters
3800 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3801 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3802 if (lipthick.gt.0.0d0) then
3803 bordliptop=(boxzsize+lipthick)/2.0
3804 bordlipbot=bordliptop-lipthick
3805 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3806 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3807 buflipbot=bordlipbot+lipbufthick
3808 bufliptop=bordliptop-lipbufthick
3809 if ((lipbufthick*2.0d0).gt.lipthick) &
3810 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3811 endif !lipthick.gt.0
3812 write(iout,*) "bordliptop=",bordliptop
3813 write(iout,*) "bordlipbot=",bordlipbot
3814 write(iout,*) "bufliptop=",bufliptop
3815 write(iout,*) "buflipbot=",buflipbot
3816 write (iout,*) "SHIELD MODE",shield_mode
3818 !C-------------------------
3819 minim=(index(controlcard,'MINIMIZE').gt.0)
3820 dccart=(index(controlcard,'CART').gt.0)
3821 overlapsc=(index(controlcard,'OVERLAP').gt.0)
3822 overlapsc=.not.overlapsc
3823 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3824 searchsc=.not.searchsc
3825 sideadd=(index(controlcard,'SIDEADD').gt.0)
3826 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3827 outpdb=(index(controlcard,'PDBOUT').gt.0)
3828 outmol2=(index(controlcard,'MOL2OUT').gt.0)
3829 pdbref=(index(controlcard,'PDBREF').gt.0)
3830 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3831 indpdb=index(controlcard,'PDBSTART')
3832 extconf=(index(controlcard,'EXTCONF').gt.0)
3833 call readi(controlcard,'IPRINT',iprint,0)
3834 call readi(controlcard,'MAXGEN',maxgen,10000)
3835 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3836 call readi(controlcard,"KDIAG",kdiag,0)
3837 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3838 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3839 write (iout,*) "RESCALE_MODE",rescale_mode
3840 split_ene=index(controlcard,'SPLIT_ENE').gt.0
3841 if (index(controlcard,'REGULAR').gt.0.0D0) then
3842 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3846 if (index(controlcard,'CHECKGRAD').gt.0) then
3848 if (index(controlcard,'CART').gt.0) then
3850 elseif (index(controlcard,'CARINT').gt.0) then
3855 elseif (index(controlcard,'THREAD').gt.0) then
3857 call readi(controlcard,'THREAD',nthread,0)
3858 if (nthread.gt.0) then
3859 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3862 write (iout,'(a)')'A number has to follow the THREAD keyword.'
3863 stop 'Error termination in Read_Control.'
3865 else if (index(controlcard,'MCMA').gt.0) then
3867 else if (index(controlcard,'MCEE').gt.0) then
3869 else if (index(controlcard,'MULTCONF').gt.0) then
3871 else if (index(controlcard,'MAP').gt.0) then
3873 call readi(controlcard,'MAP',nmap,0)
3874 else if (index(controlcard,'CSA').gt.0) then
3876 !rc else if (index(controlcard,'ZSCORE').gt.0) then
3878 !rc ZSCORE is rm from UNRES, modecalc=9 is available
3881 !fcm else if (index(controlcard,'MCMF').gt.0) then
3883 else if (index(controlcard,'SOFTREG').gt.0) then
3885 else if (index(controlcard,'CHECK_BOND').gt.0) then
3887 else if (index(controlcard,'TEST').gt.0) then
3889 else if (index(controlcard,'MD').gt.0) then
3891 else if (index(controlcard,'RE ').gt.0) then
3895 lmuca=index(controlcard,'MUCA').gt.0
3896 call readi(controlcard,'MUCADYN',mucadyn,0)
3897 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3898 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3900 write (iout,*) 'MUCADYN=',mucadyn
3901 write (iout,*) 'MUCASMOOTH=',muca_smooth
3904 iscode=index(controlcard,'ONE_LETTER')
3905 indphi=index(controlcard,'PHI')
3906 indback=index(controlcard,'BACK')
3907 iranconf=index(controlcard,'RAND_CONF')
3908 i2ndstr=index(controlcard,'USE_SEC_PRED')
3909 gradout=index(controlcard,'GRADOUT').gt.0
3910 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3911 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3912 if (me.eq.king .or. .not.out1file ) &
3913 write (iout,*) "DISTCHAINMAX",distchainmax
3915 if(me.eq.king.or..not.out1file) &
3916 write (iout,'(2a)') diagmeth(kdiag),&
3917 ' routine used to diagonalize matrices.'
3918 if (shield_mode.gt.0) then
3919 pi=4.0D0*datan(1.0D0)
3920 !C VSolvSphere the volume of solving sphere
3922 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3923 !C there will be no distinction between proline peptide group and normal peptide
3924 !C group in case of shielding parameters
3925 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3926 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3927 write (iout,*) VSolvSphere,VSolvSphere_div
3928 !C long axis of side chain
3930 ! long_r_sidechain(i)=vbldsc0(1,i)
3931 ! short_r_sidechain(i)=sigma0(i)
3932 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3938 end subroutine read_control
3939 !-----------------------------------------------------------------------------
3940 subroutine read_REMDpar
3942 ! Read REMD settings
3949 use control_data, only:out1file
3950 ! implicit real*8 (a-h,o-z)
3951 ! include 'DIMENSIONS'
3952 ! include 'COMMON.IOUNITS'
3953 ! include 'COMMON.TIME1'
3954 ! include 'COMMON.MD'
3957 !el include 'COMMON.LANGEVIN'
3959 !el include 'COMMON.LANGEVIN.lang0'
3961 ! include 'COMMON.INTERACT'
3962 ! include 'COMMON.NAMES'
3963 ! include 'COMMON.GEO'
3964 ! include 'COMMON.REMD'
3965 ! include 'COMMON.CONTROL'
3966 ! include 'COMMON.SETUP'
3967 ! character(len=80) :: ucase
3968 character(len=320) :: controlcard
3969 character(len=3200) :: controlcard1
3970 integer :: iremd_m_total
3973 ! real(kind=8) :: var,ene
3975 if(me.eq.king.or..not.out1file) &
3976 write (iout,*) "REMD setup"
3978 call card_concat(controlcard,.true.)
3979 call readi(controlcard,"NREP",nrep,3)
3980 call readi(controlcard,"NSTEX",nstex,1000)
3981 call reada(controlcard,"RETMIN",retmin,10.0d0)
3982 call reada(controlcard,"RETMAX",retmax,1000.0d0)
3983 mremdsync=(index(controlcard,'SYNC').gt.0)
3984 call readi(controlcard,"NSYN",i_sync_step,100)
3985 restart1file=(index(controlcard,'REST1FILE').gt.0)
3986 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3987 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3988 if(max_cache_traj_use.gt.max_cache_traj) &
3989 max_cache_traj_use=max_cache_traj
3990 if(me.eq.king.or..not.out1file) then
3991 !d if (traj1file) then
3992 !rc caching is in testing - NTWX is not ignored
3993 !d write (iout,*) "NTWX value is ignored"
3994 !d write (iout,*) " trajectory is stored to one file by master"
3995 !d write (iout,*) " before exchange at NSTEX intervals"
3997 write (iout,*) "NREP= ",nrep
3998 write (iout,*) "NSTEX= ",nstex
3999 write (iout,*) "SYNC= ",mremdsync
4000 write (iout,*) "NSYN= ",i_sync_step
4001 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
4004 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
4005 if (index(controlcard,'TLIST').gt.0) then
4007 call card_concat(controlcard1,.true.)
4008 read(controlcard1,*) (remd_t(i),i=1,nrep)
4009 if(me.eq.king.or..not.out1file) &
4010 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
4013 if (index(controlcard,'MLIST').gt.0) then
4015 call card_concat(controlcard1,.true.)
4016 read(controlcard1,*) (remd_m(i),i=1,nrep)
4017 if(me.eq.king.or..not.out1file) then
4018 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
4021 iremd_m_total=iremd_m_total+remd_m(i)
4023 write (iout,*) 'Total number of replicas ',iremd_m_total
4026 if(me.eq.king.or..not.out1file) &
4027 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
4029 end subroutine read_REMDpar
4030 !-----------------------------------------------------------------------------
4031 subroutine read_MDpar
4035 use control_data, only: r_cut,rlamb,out1file
4037 use geometry_data, only: pi
4039 ! implicit real*8 (a-h,o-z)
4040 ! include 'DIMENSIONS'
4041 ! include 'COMMON.IOUNITS'
4042 ! include 'COMMON.TIME1'
4043 ! include 'COMMON.MD'
4046 !el include 'COMMON.LANGEVIN'
4048 !el include 'COMMON.LANGEVIN.lang0'
4050 ! include 'COMMON.INTERACT'
4051 ! include 'COMMON.NAMES'
4052 ! include 'COMMON.GEO'
4053 ! include 'COMMON.SETUP'
4054 ! include 'COMMON.CONTROL'
4055 ! include 'COMMON.SPLITELE'
4056 ! character(len=80) :: ucase
4057 character(len=320) :: controlcard
4062 call card_concat(controlcard,.true.)
4063 call readi(controlcard,"NSTEP",n_timestep,1000000)
4064 call readi(controlcard,"NTWE",ntwe,100)
4065 call readi(controlcard,"NTWX",ntwx,1000)
4066 call reada(controlcard,"DT",d_time,1.0d-1)
4067 call reada(controlcard,"DVMAX",dvmax,2.0d1)
4068 call reada(controlcard,"DAMAX",damax,1.0d1)
4069 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
4070 call readi(controlcard,"LANG",lang,0)
4071 RESPA = index(controlcard,"RESPA") .gt. 0
4072 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
4073 ntime_split0=ntime_split
4074 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
4075 ntime_split0=ntime_split
4076 call reada(controlcard,"R_CUT",r_cut,2.0d0)
4077 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
4078 rest = index(controlcard,"REST").gt.0
4079 tbf = index(controlcard,"TBF").gt.0
4080 usampl = index(controlcard,"USAMPL").gt.0
4081 mdpdb = index(controlcard,"MDPDB").gt.0
4082 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4083 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
4084 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
4085 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
4086 if (count_reset_moment.eq.0) count_reset_moment=1000000000
4087 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
4088 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
4089 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
4090 if (count_reset_vel.eq.0) count_reset_vel=1000000000
4091 large = index(controlcard,"LARGE").gt.0
4092 print_compon = index(controlcard,"PRINT_COMPON").gt.0
4093 rattle = index(controlcard,"RATTLE").gt.0
4094 ! if performing umbrella sampling, fragments constrained are read from the fragment file
4100 if(me.eq.king.or..not.out1file) then
4102 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
4104 write (iout,'(a)') "The units are:"
4105 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
4106 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
4107 " acceleration: angstrom/(48.9 fs)**2"
4108 write (iout,'(a)') "energy: kcal/mol, temperature: K"
4110 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
4111 write (iout,'(a60,f10.5,a)') &
4112 "Initial time step of numerical integration:",d_time,&
4114 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
4116 write (iout,'(2a,i4,a)') &
4117 "A-MTS algorithm used; initial time step for fast-varying",&
4118 " short-range forces split into",ntime_split," steps."
4119 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
4120 r_cut," lambda",rlamb
4122 write (iout,'(2a,f10.5)') &
4123 "Maximum acceleration threshold to reduce the time step",&
4124 "/increase split number:",damax
4125 write (iout,'(2a,f10.5)') &
4126 "Maximum predicted energy drift to reduce the timestep",&
4127 "/increase split number:",edriftmax
4128 write (iout,'(a60,f10.5)') &
4129 "Maximum velocity threshold to reduce velocities:",dvmax
4130 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
4131 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
4132 if (rattle) write (iout,'(a60)') &
4133 "Rattle algorithm used to constrain the virtual bonds"
4137 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
4138 call reada(controlcard,"RWAT",rwat,1.4d0)
4139 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
4140 surfarea=index(controlcard,"SURFAREA").gt.0
4141 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
4142 if(me.eq.king.or..not.out1file)then
4143 write (iout,'(/a,$)') "Langevin dynamics calculation"
4145 write (iout,'(a/)') &
4146 " with direct integration of Langevin equations"
4147 else if (lang.eq.2) then
4148 write (iout,'(a/)') " with TINKER stochasic MD integrator"
4149 else if (lang.eq.3) then
4150 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
4151 else if (lang.eq.4) then
4152 write (iout,'(a/)') " in overdamped mode"
4154 write (iout,'(//a,i5)') &
4155 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4158 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4159 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4160 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4161 write (iout,'(a60,f10.5)') &
4162 "Scaling factor of the friction forces:",scal_fric
4163 if (surfarea) write (iout,'(2a,i10,a)') &
4164 "Friction coefficients will be scaled by solvent-accessible",&
4165 " surface area every",reset_fricmat," steps."
4167 ! Calculate friction coefficients and bounds of stochastic forces
4168 eta=6*pi*cPoise*etawat
4169 if(me.eq.king.or..not.out1file) &
4170 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4173 do j=1,5 !types of molecules
4174 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4175 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4177 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4178 do j=1,5 !types of molecules
4180 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4181 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4185 if(me.eq.king.or..not.out1file)then
4186 write (iout,'(/2a/)') &
4187 "Radii of site types and friction coefficients and std's of",&
4188 " stochastic forces of fully exposed sites"
4189 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4191 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4192 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4196 if(me.eq.king.or..not.out1file)then
4197 write (iout,'(a)') "Berendsen bath calculation"
4198 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4199 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4201 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4202 count_reset_moment," steps"
4204 write (iout,'(a,i10,a)') &
4205 "Velocities will be reset at random every",count_reset_vel,&
4209 if(me.eq.king.or..not.out1file) &
4210 write (iout,'(a31)') "Microcanonical mode calculation"
4212 if(me.eq.king.or..not.out1file)then
4213 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4215 write(iout,*) "MD running with constraints."
4216 write(iout,*) "Equilibration time ", eq_time, " mtus."
4217 write(iout,*) "Constraining ", nfrag," fragments."
4218 write(iout,*) "Length of each fragment, weight and q0:"
4220 write (iout,*) "Set of restraints #",iset
4222 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4223 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4225 write(iout,*) "constraints between ", npair, "fragments."
4226 write(iout,*) "constraint pairs, weights and q0:"
4228 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4229 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4231 write(iout,*) "angle constraints within ", nfrag_back,&
4232 "backbone fragments."
4233 write(iout,*) "fragment, weights:"
4235 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4236 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4237 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4240 iset=mod(kolor,nset)+1
4243 if(me.eq.king.or..not.out1file) &
4244 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4246 end subroutine read_MDpar
4247 !-----------------------------------------------------------------------------
4251 ! implicit real*8 (a-h,o-z)
4252 ! include 'DIMENSIONS'
4253 ! include 'COMMON.MAP'
4254 ! include 'COMMON.IOUNITS'
4255 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4256 character(len=80) :: mapcard !,ucase
4259 ! real(kind=8) :: var,ene
4262 read (inp,'(a)') mapcard
4263 mapcard=ucase(mapcard)
4264 if (index(mapcard,'PHI').gt.0) then
4266 else if (index(mapcard,'THE').gt.0) then
4268 else if (index(mapcard,'ALP').gt.0) then
4270 else if (index(mapcard,'OME').gt.0) then
4273 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4274 stop 'Error - illegal variable spec in MAP card.'
4276 call readi (mapcard,'RES1',res1(imap),0)
4277 call readi (mapcard,'RES2',res2(imap),0)
4278 if (res1(imap).eq.0) then
4279 res1(imap)=res2(imap)
4280 else if (res2(imap).eq.0) then
4281 res2(imap)=res1(imap)
4283 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4284 write (iout,'(a)') &
4285 'Error - illegal definition of variable group in MAP.'
4286 stop 'Error - illegal definition of variable group in MAP.'
4288 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4289 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4290 call readi(mapcard,'NSTEP',nstep(imap),0)
4291 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4292 write (iout,'(a)') &
4293 'Illegal boundary and/or step size specification in MAP.'
4294 stop 'Illegal boundary and/or step size specification in MAP.'
4298 end subroutine map_read
4299 !-----------------------------------------------------------------------------
4302 use control_data, only: vdisulf
4304 ! implicit real*8 (a-h,o-z)
4305 ! include 'DIMENSIONS'
4306 ! include 'COMMON.IOUNITS'
4307 ! include 'COMMON.GEO'
4308 ! include 'COMMON.CSA'
4309 ! include 'COMMON.BANK'
4310 ! include 'COMMON.CONTROL'
4311 ! character(len=80) :: ucase
4312 character(len=620) :: mcmcard
4314 ! integer :: ntf,ik,iw_pdb
4315 ! real(kind=8) :: var,ene
4317 call card_concat(mcmcard,.true.)
4319 call readi(mcmcard,'NCONF',nconf,50)
4320 call readi(mcmcard,'NADD',nadd,0)
4321 call readi(mcmcard,'JSTART',jstart,1)
4322 call readi(mcmcard,'JEND',jend,1)
4323 call readi(mcmcard,'NSTMAX',nstmax,500000)
4324 call readi(mcmcard,'N0',n0,1)
4325 call readi(mcmcard,'N1',n1,6)
4326 call readi(mcmcard,'N2',n2,4)
4327 call readi(mcmcard,'N3',n3,0)
4328 call readi(mcmcard,'N4',n4,0)
4329 call readi(mcmcard,'N5',n5,0)
4330 call readi(mcmcard,'N6',n6,10)
4331 call readi(mcmcard,'N7',n7,0)
4332 call readi(mcmcard,'N8',n8,0)
4333 call readi(mcmcard,'N9',n9,0)
4334 call readi(mcmcard,'N14',n14,0)
4335 call readi(mcmcard,'N15',n15,0)
4336 call readi(mcmcard,'N16',n16,0)
4337 call readi(mcmcard,'N17',n17,0)
4338 call readi(mcmcard,'N18',n18,0)
4340 vdisulf=(index(mcmcard,'DYNSS').gt.0)
4342 call readi(mcmcard,'NDIFF',ndiff,2)
4343 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4344 call readi(mcmcard,'IS1',is1,1)
4345 call readi(mcmcard,'IS2',is2,8)
4346 call readi(mcmcard,'NRAN0',nran0,4)
4347 call readi(mcmcard,'NRAN1',nran1,2)
4348 call readi(mcmcard,'IRR',irr,1)
4349 call readi(mcmcard,'NSEED',nseed,20)
4350 call readi(mcmcard,'NTOTAL',ntotal,10000)
4351 call reada(mcmcard,'CUT1',cut1,2.0d0)
4352 call reada(mcmcard,'CUT2',cut2,5.0d0)
4353 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4354 call readi(mcmcard,'ICMAX',icmax,3)
4355 call readi(mcmcard,'IRESTART',irestart,0)
4356 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
4359 call reada(mcmcard,'DELE',dele,20.0d0)
4360 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4361 call readi(mcmcard,'IREF',iref,0)
4362 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4363 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4364 call readi(mcmcard,'NCONF_IN',nconf_in,0)
4365 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4366 write (iout,*) "NCONF_IN",nconf_in
4368 end subroutine csaread
4369 !-----------------------------------------------------------------------------
4373 use control_data, only: MaxMoveType
4376 ! implicit real*8 (a-h,o-z)
4377 ! include 'DIMENSIONS'
4378 ! include 'COMMON.MCM'
4379 ! include 'COMMON.MCE'
4380 ! include 'COMMON.IOUNITS'
4381 ! character(len=80) :: ucase
4382 character(len=320) :: mcmcard
4385 ! real(kind=8) :: var,ene
4387 call card_concat(mcmcard,.true.)
4388 call readi(mcmcard,'MAXACC',maxacc,100)
4389 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4390 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4391 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4392 call readi(mcmcard,'MAXREPM',maxrepm,200)
4393 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4394 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4395 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4396 call reada(mcmcard,'E_UP',e_up,5.0D0)
4397 call reada(mcmcard,'DELTE',delte,0.1D0)
4398 call readi(mcmcard,'NSWEEP',nsweep,5)
4399 call readi(mcmcard,'NSTEPH',nsteph,0)
4400 call readi(mcmcard,'NSTEPC',nstepc,0)
4401 call reada(mcmcard,'TMIN',tmin,298.0D0)
4402 call reada(mcmcard,'TMAX',tmax,298.0D0)
4403 call readi(mcmcard,'NWINDOW',nwindow,0)
4404 call readi(mcmcard,'PRINT_MC',print_mc,0)
4405 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4406 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4407 ent_read=(index(mcmcard,'ENT_READ').gt.0)
4408 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4409 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4410 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4411 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4412 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4413 if (nwindow.gt.0) then
4414 allocate(winstart(nwindow)) !!el (maxres)
4415 allocate(winend(nwindow)) !!el
4416 allocate(winlen(nwindow)) !!el
4417 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4419 winlen(i)=winend(i)-winstart(i)+1
4422 if (tmax.lt.tmin) tmax=tmin
4423 if (tmax.eq.tmin) then
4427 if (nstepc.gt.0 .and. nsteph.gt.0) then
4428 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4429 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4431 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4432 ! Probabilities of different move types
4433 sumpro_type(0)=0.0D0
4434 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4435 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4436 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4437 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
4438 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4439 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4440 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4442 print *,'i',i,' sumprotype',sumpro_type(i)
4443 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4444 print *,'i',i,' sumprotype',sumpro_type(i)
4447 end subroutine mcmread
4448 !-----------------------------------------------------------------------------
4449 subroutine read_minim
4452 ! implicit real*8 (a-h,o-z)
4453 ! include 'DIMENSIONS'
4454 ! include 'COMMON.MINIM'
4455 ! include 'COMMON.IOUNITS'
4456 ! character(len=80) :: ucase
4457 character(len=320) :: minimcard
4459 ! integer :: ntf,ik,iw_pdb
4460 ! real(kind=8) :: var,ene
4462 call card_concat(minimcard,.true.)
4463 call readi(minimcard,'MAXMIN',maxmin,2000)
4464 call readi(minimcard,'MAXFUN',maxfun,5000)
4465 call readi(minimcard,'MINMIN',minmin,maxmin)
4466 call readi(minimcard,'MINFUN',minfun,maxmin)
4467 call reada(minimcard,'TOLF',tolf,1.0D-2)
4468 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4469 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4470 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4471 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4472 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4473 'Options in energy minimization:'
4474 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4475 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4476 'MinMin:',MinMin,' MinFun:',MinFun,&
4477 ' TolF:',TolF,' RTolF:',RTolF
4479 end subroutine read_minim
4480 !-----------------------------------------------------------------------------
4481 subroutine openunits
4483 use MD_data, only: usampl
4486 use control_data, only:out1file
4487 use control, only: getenv_loc
4488 ! implicit real*8 (a-h,o-z)
4489 ! include 'DIMENSIONS'
4492 character(len=16) :: form,nodename
4493 integer :: nodelen,ierror,npos
4495 ! include 'COMMON.SETUP'
4496 ! include 'COMMON.IOUNITS'
4497 ! include 'COMMON.MD'
4498 ! include 'COMMON.CONTROL'
4499 integer :: lenpre,lenpot,lentmp !,ilen
4501 character(len=3) :: out1file_text !,ucase
4502 character(len=3) :: ll
4505 ! integer :: ntf,ik,iw_pdb
4506 ! real(kind=8) :: var,ene
4508 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4509 call getenv_loc("PREFIX",prefix)
4511 call getenv_loc("POT",pot)
4512 call getenv_loc("DIRTMP",tmpdir)
4513 call getenv_loc("CURDIR",curdir)
4514 call getenv_loc("OUT1FILE",out1file_text)
4515 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4516 out1file_text=ucase(out1file_text)
4517 if (out1file_text(1:1).eq."Y") then
4520 out1file=fg_rank.gt.0
4525 if (lentmp.gt.0) then
4526 write (*,'(80(1h!))')
4527 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
4528 write (*,'(80(1h!))')
4529 write (*,*)"All output files will be on node /tmp directory."
4531 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4532 if (me.eq.king) then
4533 write (*,*) "The master node is ",nodename
4534 else if (fg_rank.eq.0) then
4535 write (*,*) "I am the CG slave node ",nodename
4537 write (*,*) "I am the FG slave node ",nodename
4540 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
4541 lenpre = lentmp+lenpre+1
4543 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
4544 ! Get the names and open the input files
4545 #if defined(WINIFL) || defined(WINPGI)
4546 open(1,file=pref_orig(:ilen(pref_orig))// &
4547 '.inp',status='old',readonly,shared)
4548 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4549 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4550 ! Get parameter filenames and open the parameter files.
4551 call getenv_loc('BONDPAR',bondname)
4552 open (ibond,file=bondname,status='old',readonly,shared)
4553 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4554 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
4555 call getenv_loc('THETPAR',thetname)
4556 open (ithep,file=thetname,status='old',readonly,shared)
4557 call getenv_loc('ROTPAR',rotname)
4558 open (irotam,file=rotname,status='old',readonly,shared)
4559 call getenv_loc('TORPAR',torname)
4560 open (itorp,file=torname,status='old',readonly,shared)
4561 call getenv_loc('TORDPAR',tordname)
4562 open (itordp,file=tordname,status='old',readonly,shared)
4563 call getenv_loc('FOURIER',fouriername)
4564 open (ifourier,file=fouriername,status='old',readonly,shared)
4565 call getenv_loc('ELEPAR',elename)
4566 open (ielep,file=elename,status='old',readonly,shared)
4567 call getenv_loc('SIDEPAR',sidename)
4568 open (isidep,file=sidename,status='old',readonly,shared)
4570 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4571 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
4572 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4573 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
4574 call getenv_loc('TORPAR_NUCL',torname_nucl)
4575 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
4576 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4577 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
4578 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4579 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
4582 #elif (defined CRAY) || (defined AIX)
4583 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4585 ! print *,"Processor",myrank," opened file 1"
4586 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4587 ! print *,"Processor",myrank," opened file 9"
4588 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4589 ! Get parameter filenames and open the parameter files.
4590 call getenv_loc('BONDPAR',bondname)
4591 open (ibond,file=bondname,status='old',action='read')
4592 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4593 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4595 ! print *,"Processor",myrank," opened file IBOND"
4596 call getenv_loc('THETPAR',thetname)
4597 open (ithep,file=thetname,status='old',action='read')
4598 ! print *,"Processor",myrank," opened file ITHEP"
4599 call getenv_loc('ROTPAR',rotname)
4600 open (irotam,file=rotname,status='old',action='read')
4601 ! print *,"Processor",myrank," opened file IROTAM"
4602 call getenv_loc('TORPAR',torname)
4603 open (itorp,file=torname,status='old',action='read')
4604 ! print *,"Processor",myrank," opened file ITORP"
4605 call getenv_loc('TORDPAR',tordname)
4606 open (itordp,file=tordname,status='old',action='read')
4607 ! print *,"Processor",myrank," opened file ITORDP"
4608 call getenv_loc('SCCORPAR',sccorname)
4609 open (isccor,file=sccorname,status='old',action='read')
4610 ! print *,"Processor",myrank," opened file ISCCOR"
4611 call getenv_loc('FOURIER',fouriername)
4612 open (ifourier,file=fouriername,status='old',action='read')
4613 ! print *,"Processor",myrank," opened file IFOURIER"
4614 call getenv_loc('ELEPAR',elename)
4615 open (ielep,file=elename,status='old',action='read')
4616 ! print *,"Processor",myrank," opened file IELEP"
4617 call getenv_loc('SIDEPAR',sidename)
4618 open (isidep,file=sidename,status='old',action='read')
4620 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4621 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4622 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4623 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4624 call getenv_loc('TORPAR_NUCL',torname_nucl)
4625 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4626 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4627 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4628 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4629 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4631 call getenv_loc('LIPTRANPAR',liptranname)
4632 open (iliptranpar,file=liptranname,status='old',action='read')
4633 call getenv_loc('TUBEPAR',tubename)
4634 open (itube,file=tubename,status='old',action='read')
4635 call getenv_loc('IONPAR',ionname)
4636 open (iion,file=ionname,status='old',action='read')
4638 ! print *,"Processor",myrank," opened file ISIDEP"
4639 ! print *,"Processor",myrank," opened parameter files"
4641 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
4642 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4643 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4644 ! Get parameter filenames and open the parameter files.
4645 call getenv_loc('BONDPAR',bondname)
4646 open (ibond,file=bondname,status='old')
4647 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4648 open (ibond_nucl,file=bondname_nucl,status='old')
4650 call getenv_loc('THETPAR',thetname)
4651 open (ithep,file=thetname,status='old')
4652 call getenv_loc('ROTPAR',rotname)
4653 open (irotam,file=rotname,status='old')
4654 call getenv_loc('TORPAR',torname)
4655 open (itorp,file=torname,status='old')
4656 call getenv_loc('TORDPAR',tordname)
4657 open (itordp,file=tordname,status='old')
4658 call getenv_loc('SCCORPAR',sccorname)
4659 open (isccor,file=sccorname,status='old')
4660 call getenv_loc('FOURIER',fouriername)
4661 open (ifourier,file=fouriername,status='old')
4662 call getenv_loc('ELEPAR',elename)
4663 open (ielep,file=elename,status='old')
4664 call getenv_loc('SIDEPAR',sidename)
4665 open (isidep,file=sidename,status='old')
4667 open (ithep_nucl,file=thetname_nucl,status='old')
4668 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4669 open (irotam_nucl,file=rotname_nucl,status='old')
4670 call getenv_loc('TORPAR_NUCL',torname_nucl)
4671 open (itorp_nucl,file=torname_nucl,status='old')
4672 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4673 open (itordp_nucl,file=tordname_nucl,status='old')
4674 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4675 open (isidep_nucl,file=sidename_nucl,status='old')
4677 call getenv_loc('LIPTRANPAR',liptranname)
4678 open (iliptranpar,file=liptranname,status='old')
4679 call getenv_loc('TUBEPAR',tubename)
4680 open (itube,file=tubename,status='old')
4681 call getenv_loc('IONPAR',ionname)
4682 open (iion,file=ionname,status='old')
4684 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4686 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4687 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4688 ! Get parameter filenames and open the parameter files.
4689 call getenv_loc('BONDPAR',bondname)
4690 open (ibond,file=bondname,status='old',action='read')
4691 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4692 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4693 call getenv_loc('THETPAR',thetname)
4694 open (ithep,file=thetname,status='old',action='read')
4695 call getenv_loc('ROTPAR',rotname)
4696 open (irotam,file=rotname,status='old',action='read')
4697 call getenv_loc('TORPAR',torname)
4698 open (itorp,file=torname,status='old',action='read')
4699 call getenv_loc('TORDPAR',tordname)
4700 open (itordp,file=tordname,status='old',action='read')
4701 call getenv_loc('SCCORPAR',sccorname)
4702 open (isccor,file=sccorname,status='old',action='read')
4704 call getenv_loc('THETPARPDB',thetname_pdb)
4705 print *,"thetname_pdb ",thetname_pdb
4706 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
4707 print *,ithep_pdb," opened"
4709 call getenv_loc('FOURIER',fouriername)
4710 open (ifourier,file=fouriername,status='old',readonly)
4711 call getenv_loc('ELEPAR',elename)
4712 open (ielep,file=elename,status='old',readonly)
4713 call getenv_loc('SIDEPAR',sidename)
4714 open (isidep,file=sidename,status='old',readonly)
4716 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4717 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4718 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4719 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4720 call getenv_loc('TORPAR_NUCL',torname_nucl)
4721 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4722 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4723 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4724 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4725 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4726 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
4727 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
4728 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
4729 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
4730 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
4731 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
4732 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
4733 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
4736 call getenv_loc('LIPTRANPAR',liptranname)
4737 open (iliptranpar,file=liptranname,status='old',action='read')
4738 call getenv_loc('TUBEPAR',tubename)
4739 open (itube,file=tubename,status='old',action='read')
4740 call getenv_loc('IONPAR',ionname)
4741 open (iion,file=ionname,status='old',action='read')
4744 call getenv_loc('ROTPARPDB',rotname_pdb)
4745 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
4748 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
4749 #if defined(WINIFL) || defined(WINPGI)
4750 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
4751 #elif (defined CRAY) || (defined AIX)
4752 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4754 open (iscpp_nucl,file=scpname_nucl,status='old')
4756 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4761 ! 8/9/01 In the newest version SCp interaction constants are read from a file
4762 ! Use -DOLDSCP to use hard-coded constants instead.
4764 call getenv_loc('SCPPAR',scpname)
4765 #if defined(WINIFL) || defined(WINPGI)
4766 open (iscpp,file=scpname,status='old',readonly,shared)
4767 #elif (defined CRAY) || (defined AIX)
4768 open (iscpp,file=scpname,status='old',action='read')
4770 open (iscpp,file=scpname,status='old')
4772 open (iscpp,file=scpname,status='old',action='read')
4775 call getenv_loc('PATTERN',patname)
4776 #if defined(WINIFL) || defined(WINPGI)
4777 open (icbase,file=patname,status='old',readonly,shared)
4778 #elif (defined CRAY) || (defined AIX)
4779 open (icbase,file=patname,status='old',action='read')
4781 open (icbase,file=patname,status='old')
4783 open (icbase,file=patname,status='old',action='read')
4786 ! Open output file only for CG processes
4787 ! print *,"Processor",myrank," fg_rank",fg_rank
4788 if (fg_rank.eq.0) then
4790 if (nodes.eq.1) then
4793 npos = dlog10(dfloat(nodes-1))+1
4795 if (npos.lt.3) npos=3
4796 write (liczba,'(i1)') npos
4797 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
4799 write (liczba,form) me
4800 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
4801 liczba(:ilen(liczba))
4802 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4804 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4806 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
4807 liczba(:ilen(liczba))//'.mol2'
4808 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4809 liczba(:ilen(liczba))//'.stat'
4811 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
4812 //liczba(:ilen(liczba))//'.stat')
4813 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
4816 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4817 liczba(:ilen(liczba))//'.const'
4822 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
4823 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
4824 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
4825 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
4826 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
4828 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
4830 rest2name=prefix(:ilen(prefix))//'.rst'
4832 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
4835 #if defined(AIX) || defined(PGI)
4836 if (me.eq.king .or. .not. out1file) &
4837 open(iout,file=outname,status='unknown')
4839 if (fg_rank.gt.0) then
4840 write (liczba,'(i3.3)') myrank/nfgtasks
4841 write (ll,'(bz,i3.3)') fg_rank
4842 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4847 open(igeom,file=intname,status='unknown',position='append')
4848 open(ipdb,file=pdbname,status='unknown')
4849 open(imol2,file=mol2name,status='unknown')
4850 open(istat,file=statname,status='unknown',position='append')
4852 !1out open(iout,file=outname,status='unknown')
4855 if (me.eq.king .or. .not.out1file) &
4856 open(iout,file=outname,status='unknown')
4858 if (fg_rank.gt.0) then
4859 write (liczba,'(i3.3)') myrank/nfgtasks
4860 write (ll,'(bz,i3.3)') fg_rank
4861 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4866 open(igeom,file=intname,status='unknown',access='append')
4867 open(ipdb,file=pdbname,status='unknown')
4868 open(imol2,file=mol2name,status='unknown')
4869 open(istat,file=statname,status='unknown',access='append')
4871 !1out open(iout,file=outname,status='unknown')
4874 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
4875 csa_seed=prefix(:lenpre)//'.CSA.seed'
4876 csa_history=prefix(:lenpre)//'.CSA.history'
4877 csa_bank=prefix(:lenpre)//'.CSA.bank'
4878 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
4879 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
4880 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
4881 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
4882 csa_int=prefix(:lenpre)//'.int'
4883 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
4884 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
4885 csa_in=prefix(:lenpre)//'.CSA.in'
4886 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4889 write (iout,'(80(1h-))')
4890 write (iout,'(30x,a)') "FILE ASSIGNMENT"
4891 write (iout,'(80(1h-))')
4892 write (iout,*) "Input file : ",&
4893 pref_orig(:ilen(pref_orig))//'.inp'
4894 write (iout,*) "Output file : ",&
4895 outname(:ilen(outname))
4897 write (iout,*) "Sidechain potential file : ",&
4898 sidename(:ilen(sidename))
4900 write (iout,*) "SCp potential file : ",&
4901 scpname(:ilen(scpname))
4903 write (iout,*) "Electrostatic potential file : ",&
4904 elename(:ilen(elename))
4905 write (iout,*) "Cumulant coefficient file : ",&
4906 fouriername(:ilen(fouriername))
4907 write (iout,*) "Torsional parameter file : ",&
4908 torname(:ilen(torname))
4909 write (iout,*) "Double torsional parameter file : ",&
4910 tordname(:ilen(tordname))
4911 write (iout,*) "SCCOR parameter file : ",&
4912 sccorname(:ilen(sccorname))
4913 write (iout,*) "Bond & inertia constant file : ",&
4914 bondname(:ilen(bondname))
4915 write (iout,*) "Bending parameter file : ",&
4916 thetname(:ilen(thetname))
4917 write (iout,*) "Rotamer parameter file : ",&
4918 rotname(:ilen(rotname))
4921 write (iout,*) "Thetpdb parameter file : ",&
4922 thetname_pdb(:ilen(thetname_pdb))
4925 write (iout,*) "Threading database : ",&
4926 patname(:ilen(patname))
4928 write (iout,*)" DIRTMP : ",&
4930 write (iout,'(80(1h-))')
4933 end subroutine openunits
4934 !-----------------------------------------------------------------------------
4937 use geometry_data, only: nres,dc
4939 ! implicit real*8 (a-h,o-z)
4940 ! include 'DIMENSIONS'
4941 ! include 'COMMON.CHAIN'
4942 ! include 'COMMON.IOUNITS'
4943 ! include 'COMMON.MD'
4946 ! real(kind=8) :: var,ene
4948 open(irest2,file=rest2name,status='unknown')
4949 read(irest2,*) totT,EK,potE,totE,t_bath
4952 ! AL 4/17/17: Now reading d_t(0,:) too
4954 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4957 ! AL 4/17/17: Now reading d_c(0,:) too
4959 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4962 read (irest2,*) iset
4966 end subroutine readrst
4967 !-----------------------------------------------------------------------------
4968 subroutine read_fragments
4972 use control_data, only:out1file
4975 ! implicit real*8 (a-h,o-z)
4976 ! include 'DIMENSIONS'
4980 ! include 'COMMON.SETUP'
4981 ! include 'COMMON.CHAIN'
4982 ! include 'COMMON.IOUNITS'
4983 ! include 'COMMON.MD'
4984 ! include 'COMMON.CONTROL'
4987 ! real(kind=8) :: var,ene
4989 read(inp,*) nset,nfrag,npair,nfrag_back
4991 !el from module energy
4992 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
4993 if(.not.allocated(wfrag_back)) then
4994 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4995 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4997 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4998 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
5000 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
5001 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
5004 if(me.eq.king.or..not.out1file) &
5005 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
5006 " nfrag_back",nfrag_back
5008 read(inp,*) mset(iset)
5010 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
5012 if(me.eq.king.or..not.out1file) &
5013 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
5014 ifrag(2,i,iset), qinfrag(i,iset)
5017 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
5019 if(me.eq.king.or..not.out1file) &
5020 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
5021 ipair(2,i,iset), qinpair(i,iset)
5024 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5025 wfrag_back(3,i,iset),&
5026 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5027 if(me.eq.king.or..not.out1file) &
5028 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5029 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5033 end subroutine read_fragments
5034 !-----------------------------------------------------------------------------
5036 !-----------------------------------------------------------------------------
5040 ! implicit real*8 (a-h,o-z)
5041 ! include 'DIMENSIONS'
5042 ! include 'COMMON.CSA'
5043 ! include 'COMMON.BANK'
5044 ! include 'COMMON.IOUNITS'
5046 ! integer :: ntf,ik,iw_pdb
5047 ! real(kind=8) :: var,ene
5049 open(icsa_in,file=csa_in,status="old",err=100)
5050 read(icsa_in,*) nconf
5051 read(icsa_in,*) jstart,jend
5052 read(icsa_in,*) nstmax
5053 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5054 read(icsa_in,*) nran0,nran1,irr
5055 read(icsa_in,*) nseed
5056 read(icsa_in,*) ntotal,cut1,cut2
5057 read(icsa_in,*) estop
5058 read(icsa_in,*) icmax,irestart
5059 read(icsa_in,*) ntbankm,dele,difcut
5060 read(icsa_in,*) iref,rmscut,pnccut
5061 read(icsa_in,*) ndiff
5068 end subroutine csa_read
5069 !-----------------------------------------------------------------------------
5070 subroutine initial_write
5073 ! implicit real*8 (a-h,o-z)
5074 ! include 'DIMENSIONS'
5075 ! include 'COMMON.CSA'
5076 ! include 'COMMON.BANK'
5077 ! include 'COMMON.IOUNITS'
5079 ! integer :: ntf,ik,iw_pdb
5080 ! real(kind=8) :: var,ene
5082 open(icsa_seed,file=csa_seed,status="unknown")
5083 write(icsa_seed,*) "seed"
5085 #if defined(AIX) || defined(PGI)
5086 open(icsa_history,file=csa_history,status="unknown",&
5089 open(icsa_history,file=csa_history,status="unknown",&
5092 write(icsa_history,*) nconf
5093 write(icsa_history,*) jstart,jend
5094 write(icsa_history,*) nstmax
5095 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5096 write(icsa_history,*) nran0,nran1,irr
5097 write(icsa_history,*) nseed
5098 write(icsa_history,*) ntotal,cut1,cut2
5099 write(icsa_history,*) estop
5100 write(icsa_history,*) icmax,irestart
5101 write(icsa_history,*) ntbankm,dele,difcut
5102 write(icsa_history,*) iref,rmscut,pnccut
5103 write(icsa_history,*) ndiff
5105 write(icsa_history,*)
5108 open(icsa_bank1,file=csa_bank1,status="unknown")
5109 write(icsa_bank1,*) 0
5113 end subroutine initial_write
5114 !-----------------------------------------------------------------------------
5115 subroutine restart_write
5118 ! implicit real*8 (a-h,o-z)
5119 ! include 'DIMENSIONS'
5120 ! include 'COMMON.IOUNITS'
5121 ! include 'COMMON.CSA'
5122 ! include 'COMMON.BANK'
5124 ! integer :: ntf,ik,iw_pdb
5125 ! real(kind=8) :: var,ene
5127 #if defined(AIX) || defined(PGI)
5128 open(icsa_history,file=csa_history,position="append")
5130 open(icsa_history,file=csa_history,access="append")
5132 write(icsa_history,*)
5133 write(icsa_history,*) "This is restart"
5134 write(icsa_history,*)
5135 write(icsa_history,*) nconf
5136 write(icsa_history,*) jstart,jend
5137 write(icsa_history,*) nstmax
5138 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5139 write(icsa_history,*) nran0,nran1,irr
5140 write(icsa_history,*) nseed
5141 write(icsa_history,*) ntotal,cut1,cut2
5142 write(icsa_history,*) estop
5143 write(icsa_history,*) icmax,irestart
5144 write(icsa_history,*) ntbankm,dele,difcut
5145 write(icsa_history,*) iref,rmscut,pnccut
5146 write(icsa_history,*) ndiff
5147 write(icsa_history,*)
5148 write(icsa_history,*) "irestart is: ", irestart
5150 write(icsa_history,*)
5154 end subroutine restart_write
5155 !-----------------------------------------------------------------------------
5157 !-----------------------------------------------------------------------------
5158 subroutine write_pdb(npdb,titelloc,ee)
5160 ! implicit real*8 (a-h,o-z)
5161 ! include 'DIMENSIONS'
5162 ! include 'COMMON.IOUNITS'
5163 character(len=50) :: titelloc1
5164 character*(*) :: titelloc
5165 character(len=3) :: zahl
5166 character(len=5) :: liczba5
5168 integer :: npdb !,ilen
5172 ! real(kind=8) :: var,ene
5176 if (npdb.lt.1000) then
5177 call numstr(npdb,zahl)
5178 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5180 if (npdb.lt.10000) then
5181 write(liczba5,'(i1,i4)') 0,npdb
5183 write(liczba5,'(i5)') npdb
5185 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5187 call pdbout(ee,titelloc1,ipdb)
5190 end subroutine write_pdb
5191 !-----------------------------------------------------------------------------
5193 !-----------------------------------------------------------------------------
5194 subroutine write_thread_summary
5195 ! Thread the sequence through a database of known structures
5196 use control_data, only: refstr
5198 use energy_data, only: n_ene_comp
5200 ! implicit real*8 (a-h,o-z)
5201 ! include 'DIMENSIONS'
5203 use MPI_data !include 'COMMON.INFO'
5206 ! include 'COMMON.CONTROL'
5207 ! include 'COMMON.CHAIN'
5208 ! include 'COMMON.DBASE'
5209 ! include 'COMMON.INTERACT'
5210 ! include 'COMMON.VAR'
5211 ! include 'COMMON.THREAD'
5212 ! include 'COMMON.FFIELD'
5213 ! include 'COMMON.SBRIDGE'
5214 ! include 'COMMON.HEADER'
5215 ! include 'COMMON.NAMES'
5216 ! include 'COMMON.IOUNITS'
5217 ! include 'COMMON.TIME1'
5219 integer,dimension(maxthread) :: ip
5220 real(kind=8),dimension(0:n_ene) :: energia
5222 integer :: i,j,ii,jj,ipj,ik,kk,ist
5223 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5225 write (iout,'(30x,a/)') &
5226 ' *********** Summary threading statistics ************'
5227 write (iout,'(a)') 'Initial energies:'
5228 write (iout,'(a4,2x,a12,14a14,3a8)') &
5229 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5230 'RMSnat','NatCONT','NNCONT','RMS'
5231 ! Energy sort patterns
5236 enet=ener(n_ene-1,ip(i))
5239 if (ener(n_ene-1,ip(j)).lt.enet) then
5241 enet=ener(n_ene-1,ip(j))
5253 ist=nres_base(2,ii)+ipatt(2,i)
5255 energia(i)=ener0(kk,i)
5257 etot=ener0(n_ene_comp+1,i)
5258 rmsnat=ener0(n_ene_comp+2,i)
5259 rms=ener0(n_ene_comp+3,i)
5260 frac=ener0(n_ene_comp+4,i)
5261 frac_nn=ener0(n_ene_comp+5,i)
5264 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5265 i,str_nam(ii),ist+1,&
5266 (energia(print_order(kk)),kk=1,nprint_ene),&
5267 etot,rmsnat,frac,frac_nn,rms
5269 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5270 i,str_nam(ii),ist+1,&
5271 (energia(print_order(kk)),kk=1,nprint_ene),etot
5274 write (iout,'(//a)') 'Final energies:'
5275 write (iout,'(a4,2x,a12,17a14,3a8)') &
5276 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5277 'RMSnat','NatCONT','NNCONT','RMS'
5281 ist=nres_base(2,ii)+ipatt(2,i)
5283 energia(kk)=ener(kk,ik)
5285 etot=ener(n_ene_comp+1,i)
5286 rmsnat=ener(n_ene_comp+2,i)
5287 rms=ener(n_ene_comp+3,i)
5288 frac=ener(n_ene_comp+4,i)
5289 frac_nn=ener(n_ene_comp+5,i)
5290 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5291 i,str_nam(ii),ist+1,&
5292 (energia(print_order(kk)),kk=1,nprint_ene),&
5293 etot,rmsnat,frac,frac_nn,rms
5295 write (iout,'(/a/)') 'IEXAM array:'
5296 write (iout,'(i5)') nexcl
5298 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5300 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5301 'Max. time for threading step ',max_time_for_thread,&
5302 'Average time for threading step: ',ave_time_for_thread
5304 end subroutine write_thread_summary
5305 !-----------------------------------------------------------------------------
5306 subroutine write_stat_thread(ithread,ipattern,ist)
5308 use energy_data, only: n_ene_comp
5310 ! implicit real*8 (a-h,o-z)
5311 ! include "DIMENSIONS"
5312 ! include "COMMON.CONTROL"
5313 ! include "COMMON.IOUNITS"
5314 ! include "COMMON.THREAD"
5315 ! include "COMMON.FFIELD"
5316 ! include "COMMON.DBASE"
5317 ! include "COMMON.NAMES"
5318 real(kind=8),dimension(0:n_ene) :: energia
5320 integer :: ithread,ipattern,ist,i
5321 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5323 #if defined(AIX) || defined(PGI)
5324 open(istat,file=statname,position='append')
5326 open(istat,file=statname,access='append')
5329 energia(i)=ener(i,ithread)
5331 etot=ener(n_ene_comp+1,ithread)
5332 rmsnat=ener(n_ene_comp+2,ithread)
5333 rms=ener(n_ene_comp+3,ithread)
5334 frac=ener(n_ene_comp+4,ithread)
5335 frac_nn=ener(n_ene_comp+5,ithread)
5336 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5337 ithread,str_nam(ipattern),ist+1,&
5338 (energia(print_order(i)),i=1,nprint_ene),&
5339 etot,rmsnat,frac,frac_nn,rms
5342 end subroutine write_stat_thread
5343 !-----------------------------------------------------------------------------
5345 !-----------------------------------------------------------------------------
5346 end module io_config